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Chapter 1 
Introduction 


1.1 General Comments About This Thesis 


Each chapter of this thesis contains work that has already being published as a conference or journal 
paper, and/or is being considered for journal publication. Because of this, each chapter has an introduction 
that motivates the particular aspects of the work relevant to that chapter. Therefore, in this introductory 
chapter, a general but brief overview that is relevant to the thesis topic is presented. 


1.2 Introduction 


Joining dissimilar materials is an important and contemporary subject receiving great attention from the 
aerospace engineering sector due to the large number of structural components that are being produced with 
a variety of materials, including metals, polymers, fiber reinforced composites (FRCs) and structural foams 
due to increasing demands on weight reduction. Traditional methods of joining, which have been largely 
based on bolted joint technology, are not very suitable for joining FRCs to each other or FRCs to other 
materials [1]. 

As bonded joints increase in popularity and use, the demand for modeling techniques increases also. In 
the past, analytical models have been favored as the preferred method of predicting stresses and strength, but 
finite element (FE) methods have emerged as the new standard in preliminary design due to necessity of 
analyzing and designing components that contain multiple joints where analytical techniques become 
intractable. FE based methods have been proven to be extremely powerful, but the small scale of the 
adhesive thickness when compared to the dimensions of the surrounding structure has kept joint FE analysis 
largely out of global vehicle models. A fine mesh is needed to correctly model the adhesive layer producing 
an incompatibility in simultaneously analyzing the joint stresses accurately in conjunction with a very coarse 
model of an entire vehicle. Therefore, the actual design and sizing of joints is often put off until a later time, 
when small sub-models are used to look into the details of a vehicle. 

To further complicate things, the eccentricity of the load path in most aerospace structural joints (such 
as the popular single lap joint [2] featured in Figure 2-5) causes significant bending when axially loaded. 
This bending-introduced rotation makes nonlinear geometric effects significant, even in early stages of 
loading. Furthermore, most modern advanced polymeric adhesives show considerable nonlinear material 
behavior, which causes the joint to remain intact even when the yield stress of the adhesive is reached. 
Therefore, it is imperative that that these factors be included in models in order to fully utilize the superior 
capability of bonded joints. 

Motivated by these reasons, the overarching objective of this thesis was to develop a single finite 
element that can capture accurately the stress and strain states of a bonded joint while still facilitating its 
merger with the surrounding structure without having an incompatibility in finite element mesh densities. 
The bonded joint element (simply referred to as "joint element") uses analytical structural models to find 
shape functions for a joint, allowing the joint region to be modeled with one element. That way, with a 
single element modeling the joint, it can be inserted into a larger finite element model efficiently. 
Furthermore, large rotations, nonlinear material properties, and crack growth capabilities are included and 
analysis guidelines are developed to keep the number of elements required to a minimum. 


1.3 Thesis Organization 


The chapters are based on previously published papers, and can be read as stand-alone pieces of work. 
For this reason, the introduction and motivation for each chapter might seem repetitive. However, in the 
actual body of the chapter, many aspects have been removed from their previous state to shorten up the 
work and avoid repetition, especially in the formulation. The notation has been largely unified, with one 
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difference being that in the earlier chapters, there is only one coordinate system. In the final chapter before 
the conclusion, the co-rotational formulation necessitates a local and global coordinate system. This local 
coordinate system correlates with the coordinate system in the earlier two chapters. 


1.4 Significant Contributions 


This section provides an executive summary of the significant contributions of this dissertation to the 
scientific community. 


A linear elastic joint finite element used to find the stresses, strains, and displacements in an 
adhesively boned joint. 

Method of modeling complex, modern joints with tapers, ply drops and angle changes using 
discrete joint element “building blocks.” 

A joint element that can model the performance of a joint with a functionally graded adhesive with 
just one element. 

Demonstrated the stress reduction, sensitivity, and universal applicability of functionally graded 
adhesives. 

A technique for applying a co-rotational formulation to layered beam models to capture large 
rotation effects. 

A technique of internally re-meshing an element to introduce and grow a crack within the element. 
A technique of adapting the shape functions in an element as the loading increases to account for 
material softening and improve elemental mesh convergence. 

A method of using bulk adhesive tensile data in an adhesive spring model with uncoupled shear 
and peel responses. 
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Chapter 2 
Linear Elastic Modeling of Adhesively Bonded Joints Using an 
Enhanced Joint Finite Element 


The design and sizing of adhesively bonded joints has always been a major bottleneck in the design of 
composite vehicles. Dense meshes are required to capture the full behavior of a joint, but these dense 
meshes are impractical in vehicle-scale models where a course mesh is more desirable to make quick 
assessments and comparisons of different joint geometries. Analytical models are often helpful in sizing, but 
difficulties arise in coupling these models with full-vehicle finite element (FE) models. Therefore, a reduced 
order joint finite element was created that can be used within structural FE models to make quick 
assessments of bonded composite joints. The shape functions of the joint finite element were found by 
solving the governing equations for a joint where the adherends were modeled as beams connected by 
adhesive modeled using various assumptions. By analytically determining the shape functions of the joint 
element, the complex joint behavior can be captured with very few elements. Analyses of joint stresses for 
different joints using the enhanced joint finite elements were found to agree well with analyses using 
standard, 2-D plane stress elements. 


2.1 Introduction 


With the increasing demand for composites in lightweight aerospace structures, adhesively bonded 
joints are becoming more critical than ever. Bolts and rivets introduce holes which cause significant stress 
concentrations and premature failure in composite materials, while adhesives spread the load more evenly 
over the composite while facilitating a lighter overall structure. 

Traditionally, analytical models have been used to assess the performance of joints [2-7]. However, FE 
modeling has emerged as a popular and robust method for structural analysis. In order to properly predict 
the stresses in a joint using FE analysis, a dense mesh must be used, particularly near the adhesive ends. 
While such a model is ideal for detailed analysis, it can be crippling when it comes to joint design and sizing 
[8,9]. Furthermore, such a model does not couple easily with coarse, vehicle-scale models used for sizing of 
the vehicle. Therefore, a need exists to develop predictive tools for bonded joints that can be seamlessly 
coupled with large scale structural analyses without adding major computational demands. Such tools can 
be used to make quick mesh-independent assessments of bonded composite joints. Furthermore, they fit in 
into the computational hierarchy of virtual testing of aircraft structures [10], an area that is getting increased 
attention in the aerospace industry with the aim of lowering design cycle and certification costs. 

A solution to this problem involves merging analytical models with finite elements. Simplified 
structural models can be used to obtain shape functions that are exact for the assumptions of the model. 
These shape functions can be used to formulate stiffness matrix for the problem at hand. As long as the 
assumptions remain valid, such an element would give the exact solution regardless of the number of 
elements used. 

This method has been used to calculate an stiffness matrix for different beam on elastic foundation 
problems [11,12]. More recently, Waas and Gustafson [13] have created an element to capture the behavior 
of a double overlap joint subjected to mechanical and thermal loads. The current authors have continued this 
effort by creating an stiffness matrix for a single lap joint [14]. 

However, joints in application rarely resemble the simple joints that can be solved easily to obtain 
enhanced shape functions. Adherend tapers and steps do not conform well to a constant thickness adherend 
model. The current study lays out the formulation for the stiffness matrix for a joint with an arbitrary 
number of adherends. Additionally, an approach is introduced to model realistic, complex joints using 
simple joint elements as building blocks. This approach is demonstrated using two different joint types, and 
compared with dense mesh 2-D element models to show: 1) how well the model predicts stress near an 
adherend step and 2) how well a smooth taper can be modeled using stepped joint elements. 
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Figure 2-1. Geometric parameters for overlap region of an adhesively bonded joint 
with multiple bonded layers and the width of the joint in the y-direction is b. 


2.2 Analytical Formulation 


The enhanced joint finite element uses an analytical formulation to obtain the stiffness matrix for N 
number of adherends joined by N-1 adhesive layers. The geometric parameters around the ith adhesive layer 
are shown in Figure 2-1. The adhesives and adherends were assumed to be linearly elastic, but not 
necessarily isotropic. It is also assumed that the adherend stiffness is much greater than that of the adhesive, 
which is typical of standard aerospace applications. The strains and rotations are considered small. The 
subscript i refers to values associated with adherend i, and ai refers to adhesive layer i. Each adhesive and 
adherend layer has its own z-coordinate starting at the centerline of the layer, marked as z,; and z; 
respectively. The variables u,(x) and w(x) refer to the x and z-direction centerline displacements of 
adherend i. 

The first section below shows different models considered for the adherends, namely an Euler-Bernoulli 
and a Timoshenko model. The next section lays out various models utilized to represent the adhesive layers. 
Using various combinations of adherend and adhesive models, a system of governing equations can be 
formulated. Regardless of the models chosen, this system of governing equations can be solved following 
the method outlined in the third section, and the stiffness matrix can be obtained. Subsequent discussion will 
consider the different models and why or under what circumstances their use is favorable. 


2.2.1 Adherend Models 


The models for the adherends utilize assumptions for the stress, strain, and displacement in order to 
simplify the problem to make the possibility of obtaining an analytical solution within reach. For all models, 
the adherends are assumed to be linearly elastic, and layered transversely isotropic. It is assumed that the 
beam is very long in the y-direction, which puts the adherends in a state of cylindrical bending. The first 
model assumes the adherends to behave as an Euler-Bernoulli beam (with the assumption of plane strain 
rather than plane stress), while the second model includes shear deformations. Rather than presenting a 
detailed derivation for the three beam theories, three vectors/matrices needed to use these models will be 


defined for each model: D, , G,, and u,. The matrix D,;‘ is defined with the relation 
ot = Dz, 2-1 


where 
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The normal stress and strain are in the x-direction of the kth layer of the ith adherend, while the shear stress 
and strain are in the xz plane. Similarly, the matrix G, relates the strain to the centerline displacements with 
the equation 


€,=Guu, ee 
where U, , the centerline displacements of the ith adherend, will be defined for each model. Now D* ; G, : 
and u, can be defined for each model, enabling a unified notation and simplifying a change between the 
beam theories. 

2.2.1.1 Euler-Bernoulli Model 


The Euler-Bernoulli beam model assumes that the axial stress due to extension and bending is the most 


dominant stress, and that all others can be ignored. The stress in layer k is related to the strain by D 7 , which 
is given by: 


—k 
pi=|@u 9 ae 
‘ 0 0 


where Q,, is the 1,1 component of the transformed lamina stiffness matrix, Q_, as defined by Classical 


Lamination Theory [11,12]. The adherend strain is related to the adherend centerline displacements by G, , 
given as: 


1 


BG: 20.0. <O- °0 


[O° 1 0 0 =z; 0 2-6 


Furthermore, the centerline displacement vector is defined as 


2-7 
u,=[u(x) U(X), Wi) WOOD. Wm WD ux | 


where ,x denotes the derivative with respect to x. This form of defining the centerline displacements might 
not be conventional, but it is used to lead into our solution strategy of the governing equations. Using state 
variables with higher order derivatives as is done here allows the governing equations to be reduced to a 
series of first order differential equations. 
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2.2.1.2 Shear Deformation Model 
The shear deformation model assumes that plane sections remain plane, but not necessarily 


perpendicular to the centerline. The stress in layer k is related to the strain by D é , which is given by: 


xk 2-8 
= Qu 0 
—k 
0 kKQ¢6. 


D+ 


1 


where the variable k is a correction factor used to offset the error caused by the simplifying assumption that 
the shear stress is uniform throughout a cross section of the beam. For a rectangular cross-section, the value 
of 5/6 is most commonly used [15]. The adherend strain is related to the adherend centerline displacements 


by G, , given as: 


fo100 02 


i 


i= 2-9 
0001-1 0 
Furthermore, a new variable is introduced to define the centerline displacements. The rotation due to 
bending only is defined by W;(X) , making the centerline displacements: 
2-10 


T 


u; =[u;(x) U;(X),. w;(X) w, (x) , W(x) v(x), | 


2.2.2 Adhesive Models 


The adhesive models considered differ in the amount of detail included in the stress-strain relationship 
and the strain-displacement relationship. However, all models shown here will start out with the assumption 
that the displacements vary linearly in the z-direction. This allows the formulation to be strictly in terms of 
functions of x, as with the adherends. The adhesive is assumed to be perfectly bonded to the adherends. The 


difference in the models will be defined by two matrices: D,; , and G,,;. Similar to the adherends, the 


ai? 


matrix D_, relates the adhesive stresses to strains through 


= 2-11 
6, =D, 


ai” ai 


where 


if 2-12 


and 


fai = [Ec Vai Exai l' ; ou 


The normal stress and strain, O,,, and € 


ai? 


are the stresses and strain in the z-direction, also commonly 


referred to as the peel stress and strain. The shear stress and strain, T,, and 7,; , are in the xz plane, while the 
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axial stress and strain, O|,,, and €,,,; occur in the x-direction. Although this might not be the conventional 


i 
way of ordering stresses and strains in elasticity, this is the chosen scheme that adheres closer to conventions 
when dealing with adhesive layers. The peel and shear stresses are the components of most concern, and the 
axial stress and strain are often ignored [3,4,6,5]. Therefore, it is placed at the end to allow for easy deletion 

if desired. 


The matrix G_,; relates the strain of the adhesive to the adhesive displacements with the equation 


G 2-14 


where the adhesive displacements and derivatives, U,,, are defined as 


ai? 


U(X, Z,i) as 


UG (x, Lai ) x 


Wai (x, Z,i) 


ai 


Wai (x, Lai ys 


The adhesive displacements are related to the adherend centerline displacements of the adherends above and 
below the adhesive layer through the equation 


_. fu, 2-16 
U, - H,; 


where 
Cc 0 0 000cC, 0 0 000 oe 
Be | HOP KE h OP > ONO: ON GOA. Oo oO, 0D 
““]0 0 C€ 000 0 0c 000 
0 0 0C 00 0 0 0 C, 00 
and 
; 2-18 
Gee lit, “Ste ih 


This relation was found by assuming that the displacements in the adhesive layer vary linearly in the x- 
direction, and that the adhesive and adherends are bonded perfectly. Using the above relations, the adhesive 
strain can be related to the adherend centerline displacements by 


u, 2-19 
© ai = G,, 
Ui, 


where 
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G.=-GH. 2-20 


al al al 


G ,; will be defined for each adhesive model in the following sections. 


2.2.2.1 Model 1: Winkler Foundation 


In this model, it is assumed that the adhesive is a bed of uncoupled linear shear and normal springs. The 
adhesive layer is assumed to be so thin that the stress in the adhesive layer is independent of the z- 
coordinate. The stress is related to the strain by the matrix 


E. 0 O 2-21 
D,; = 0 G,; 0 
0 0 O 


where E,,; and G,, are the normal and shear moduli of the ith adhesive layer, respectively. As can be seen, 


the axial stress is neglected which is a common assumption for adhesive layer models [4—6,16]. The normal 
strain is defined in the standard small strain elasticity manner, but the shear strain is simplified to 


_ OU, (% Z,;) 2-22 


Vai oz 


ai 


which makes it a true Winkler foundation, or bed of uncoupled linear shear springs. From this 
simplification, the matrix relating the adhesive strains with the adherend centerline displacements, G,,, is 


given as 
0.0.2 0. 00. 0 OS oO 00 2-23 
G,=|7 00 ~~ 00700 #* 0 0 


000 0 00 0 0 0 0 0 0 


This model is the simplest model, and will be utilized heavily in the study. The main reason for this is not 
necessarily its simplicity, but lies in the fact that it resembles a Cohesive Zone Model (CZM), with peel and 
shear components uncoupled. Although a significant simplification, the model is “tuned” using experimental 
tests to characterize the peel and shear components (Mode I and Mode IJ) and is capable of predicting 
mixed-mode failure events. 


2.2.2.2 Model 2: Coupled Springs 


The second model resembles the first model, except in the strain-displacement relation. The full shear 
strain equation is used: 


2-24 
a Ou, ;(X; Z,;) ‘ Ow, (X; Zai) 
OZ; OX,,; 


al 
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This is the equivalent of having a non-Winkler foundation as the adhesive, where the normal springs are 
“tied together” and add to the shear rigidity. This model makes the adhesive shear strain/stress a function of 
Z, as can be seen in 


00+ 0 00002 0 0 0 2-25 


Nai Nai 


Gi) 2) 00. SEO Oa. 0 SG, '0: 0 


a Nai 2Nai Nai 2Nai 


OF 0 0 Oe ORB eB OO 


2.2.2.3 Model 3: Plane Stress 


Model 3 uses the full stress-strain and strain/displacement relations for a body in plane stress. This 
model is appropriate when the depth of the joint (b) is much smaller than the length of the adherends. The 
stress for the plane stress problem is related to the strain by the matrix 


Ur wee 2-26 


where V_, is the Poisson’s Ratio of the ith adhesive layer, and 


Ee 2-27 
aang 


The adhesive strain is related to the adherend centerline displacements by G,, given as: 


00 - 0 0000-2 0 0 0 2-28 
G,=|t 0 0 ~+C 0 0S 0 0 HH+C, 0 0 
0 CG 0 0 {C 0 0 Cc, 0 0 SEG. 0 


2.2.2.4 Model 4: Plane Strain 


Model 4 is similar to model 3, except that the stress-strain relation is for a body in plane strain. 
Although this model was not used in the results, it is important for modeling the joining of wide panels. The 
stress for the plane strain problem is related to the strain by the matrix 


Cy (l—- Va ViiCa aa 
D, = 0 G,; 0 
Vai C, 0 C, (1 Vai ) 


where V_, is the Poisson’s Ratio of the ith adhesive layer, and 
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. 2-30 
ai 


: (1+V,, )(1-2v4) 


2.2.3 General Formulation 


Regardless of the adherend and adhesive models chosen, the strain energy of the joint, Ujoin, can be 
written as: 


2-31 


1 N M N-1 
U cint = a bs 2 Jia of ¢dA om bs 6," 24dA i 


where At is the area (in the yz plane) of the kth layer of the ith adherend, A,, is the area of the ith adhesive 


layer, and M is the number of transversely isotropic layers of the ith adherend, and all other variables are 
previously defined. Using the principle of stationarity of potential energy, 2N fully coupled governing 
equilibrium differential equations are obtained from the energy expression. Of the 2N governing equations, 
N equations correspond to the axial equilibrium, while N equations correspond to the transverse equilibrium. 
The axial displacement equilibrium equations contain second order derivatives, while the transverse 
displacement equations have fourth order derivatives. The order of these equations can be reduced and 
assembled into a system of first order constant coefficient homogeneous ordinary differential equations of 
the form 


2-32 


where 


u=[u," K u,” K uy’ | 2-33 


or a vector containing the centerline displacements of all N adherends in the overlap region. 

Inspecting the matrix A can be helpful in determining the nature of the solution and determining the 
solution method. There are 6N eigenvalues of A: N real eigenvalues, 2N complex eigenvalues, and 3N 
repeating eigenvalues. Therefore, the solution is made up of N exponential terms, 2N exponential terms 
multiplied by a sine or cosine, and the 3N repeating eigenvalues correspond to a third order polynomial 
found in a standard beam solution. Such a complex solution shows that merely employing standard beam 
shape functions to the joint problem would be inadequate in capturing the nature of the whole solution. 

The system in can be solved using various methods, but calculating the matrix exponential was the 
chosen method because numerical boundary conditions are not required to obtain a solution. The solution of 


the system can be written in terms of the matrix exponential, e* anda vector of unknown constants, C , as 


u=e"C_ ag 
The matrix exponential can be expressed as the infinite series [17] 
M co yk : 2-35 
ev = Zs a" 
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In order to obtain faster convergence, a method of scaling and squaring [18] was employed, and the 
series was calculated up to a value of k which yields an acceptable error, € . The error can be defined many 
ways, but the current study defined the error as the difference between the 1-norms of e™ for k-1 and k. The 
value of the acceptable error was set at € = 0.0001. 

The next step was to solve for the vector of constants, C , using the boundary conditions. This is where 
the analytical formulation is discretized, and the displacements are obtained in terms of the nodal 
displacements as defined in Figure 2-2. For adherend i, the boundary conditions on the left side of the joint 
(x=0) can be expressed in the following equation: 


q' =b.u,(0) ™ 


where u,(0) is u, evaluated at x=0, q' is a vector containing the prescribed nodal degrees of freedom of 
adherend i at x=0 (or the left side), and 


0 0 2-37 
0 0 
0 O° 


Equation 2-36 for all N adherends can be assembled together, and a relation between the nodal degrees of 
freedom at x=0 and the vector of constants can be found using: 


2-38 


q = Bye*'C, 


Figure 2-2. Boundary conditions for adherend i: prescribed nodal displacements and rotations at x=0 and x=. 
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where 


u 2-39 


and 


b 2-40 


where the subscript N denotes the number of matrices on the diagonal. After performing the same operations 
at x=l, all of the boundary conditions can be combined in the form 


l 2-41 
ZC, = ("| =q. 
q 


AO 2-42 
B e 
Z= 
2N eo! 


Using this relation, one can obtain an expression for the vector of unknown constants: 


where 


sheath 2-43 
C,=z q 


This relation can be inserted into Equation 2-34 to get the adherend centerline displacements in terms of the 
nodal degrees of freedom, 


u=Nq 2-44 


where the shape functions, N are defined as: 


: 2-45 
N=e"z" 


Next, Equation 2-31 is rewritten in terms of the centerline displacements using the stress-strain and strain- 
displacement relations and put into matrix form: 


NASA/CR—2012-217606 12 


Single Lap 


Double Lap 


Figure 2-3. The building block approach facilitates modeling complex joints with 
simple joint element building blocks. 


l 2-46 
_ T . 
U joint = i u Uudx 


The usefulness of this form is that the integral has been performed analytically over the area since the 
variation in the z-direction and y-direction is known, and numerical integration must only be performed in 
the x-direction. Inserting Equation 2-24 into Equation 2-46 and minimizing the energy yields the stiffness 
matrix of the joint: 


Cs 2-47 
K jin = |, N'UNax 


This formulation gives the stiffness matrix for a simple region of constant thickness adherend overlap. 
However, many joints in application contain complicated geometries, including ply steps and tapers. To use 
the simple joint element for efficient modeling of complex joints, a building block approach was 
implemented. This approach involves combining simple, constant-thickness joint sections to create 
complicated joints with very few elements. This concept is illustrated in Figure 2-3, where single, double 
and triple adherend joint building blocks are combined to make complicated joints such as a pi joint, tapered 
single lap joint, and spliced sandwich joint. 


Adherends joined together at the adherend centerlines are related to each other with the equation 


2-48 
dit = Tir | i=1..3 


where the subscript L is for the adherend on the left, R is for the adherend on the right, and the numerical 
subscripts 1, 2, and 3 refer to axial, transverse, and rotational degrees of freedom respectively. To model a 
ply step or taper, the transverse and rotational degrees of freedom are equal, but the axial degree of freedom 
of the left adherend is related to that of the right adherend through the following equation: 


2-49 
Git = Tir ~ Coppser Isr 


where t_,.,, is the vertical (z-direction) offset distance between the two nodes (Figure 2-4). A discussion of 


offset 
the accuracy of this approximation will be discussed subsequently. 
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Figure 2-4. Joints connected together with an offset between centerlines are related using Eq. (21). 


2.3 Validation and Results 


2.3.1 Adhesive Model Comparison 


Two studies were done to validate the joint element and discover the simplest adhesive model required 
to accurately predict the shear and peel stress in a single lap joint. First, the adhesive models were compared 
with each other over a broad range of joint geometric parameters to show which models have an impact on 
the predicted stress for different parameter values. Second, a 2-D solution, based on the finite element 
method (FEM) was generated for four parameter cases and the predicted adhesive stresses were compared 
with the three models. This was done to illustrate the accuracy of the joint element for different geometric 
parameters. 

The three adhesive models were compared over a range of parameters to determine the difference in 
peak shear and peel stress along the adhesive centerline (z,;=0) predicted by the models as a function of the 
parameters. This is useful to show when assumptions about the adhesive stresses are valid and when one 
should be careful using them. Figure 2-5 shows the dimensions and material properties of the single lap joint 
considered. The adherends were aluminum, with a Young’s modulus of 70 GPa (EF) and Poisson’s ratio of 
0.33 (v; and v2). They were both 5 mm thick (t), 2 mm wide (b), and extended 5000 mm long past the joint 
overlap (ladherend). The adhesive was FM300 and had a Young’s modulus of 2.17 GPa (E,) with a shear 
modulus of 0.89 GPa (G,). The adhesive thickness (7) was varied from 0.005 to 5 mm, and the overlap 
length (1) was varied from 50 to 5000 mm. Although these parameters are not necessarily typical for a joint, 
it was necessary to test a wide range of parameters to demonstrate the working range of the joint element 
and the limits of the assumptions. 

On either side of the joint element, 50 beam elements were used to model the adherends outside of the 
overlap region. The left end was clamped, or restrained from displacement. The right end was extended in 
the x-direction by 10 mm (A) and restrained from rotating and from displacing in the z-direction as shown in 
Figure 2-5. 

Along with comparing the three models to each other, they were also compared with the results from a 
2-D FEM created using the commercial package Abaqus [19] for four different geometric parameter cases. 
The same geometric parameters and boundary conditions shown in Figure 2-5 were used for the cases, and 
Table 2-1 shows the values of t/] and 7/t used for each case. Cases were chosen at four extreme corners of 
the parameters tested for the comparison study. 
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Figure 2-5. Boundary conditions, loading scenario and geometric parameters for the single lap joint. 
Material properties include: E,;=E,= E, E,, and G,. 


Table 2-1. Parameters of in-depth study cases. 


t/l nit 
Case 1 0.001 1 
Case 2 0.001 0.1 
Case 3 0.1 1 
Case 4 0.1 0.1 


The finite element models were constructed using 70,000 - 100,000 2-D plane stress quadrilateral 
elements, with an element bias towards the corners of the adhesive. One issue that deserves mention is that 
the joint element models do not fulfill the traction free (a, and t, ) boundary condition at the free edges of 
the adhesive boundary. Therefore, the FEMs and the joint element models are not expected to predict 
similar stresses at the free edges of the adhesive. Moreover, the inside corners of the adhesive in the FEM 
causes a Stress singularity, making the model mesh dependant in the corner singularity region. In 
application, spew fillets are present, which eliminate this stress concentration and cause the edges of the 
overlap to not have zero stress [20]. Therefore it is of no major concern that the joint element does not 
reflect the stresses predicted by the FEA model at the ends of the adhesive. Since the model is meant to 
serve as an initial vehicle-scale model element for initial sizing and not necessarily for detailed analysis, the 
goal of the validation is to show that the overall behavior of the joint is reflected by the joint element. 

The plots comparing the relative difference between the maximum centerline peel and shear stresses 
predicted by Models 3 and 2 for different 7/t and t/I values can be found in Figure 2-6. The relative 
difference between the maximum peel stress predicted by Models 3 and 2 is between 8% and 5% for the 
whole range of adhesive thicknesses and overlap lengths considered. The difference is greatest for the very 
thin adhesive layers, and least for the thickest adhesive layers. Additionally, the difference between the 
maximum adhesive centerline shear stress predicted by the two models is below 1% for all of the parameters 
considered in this study. Considering the uncertainty of the stress state at the corner of the adhesive due to 
the previously mentioned issues related to stress concentrations and violation of the traction free condition, 
the difference between these models is very small. If one is only concerned with the peel and shear stress in 
the adhesive layer, it would be advantageous to use Model 2 over Model 3 due to its increased simplicity. 
However, Tsai and Morton [21] note that the magnitude of o,, near the end of the overlap is often 
comparable to the magnitude of the peel and shear stresses. Therefore, it should be kept in mind that Model 
3 is the only one of the three models to yield the extensional stress, o4,.. 

Figure 2-7 contains a comparison of the relative difference between the maximum centerline peel and 
shear stresses predicted by Models 2 and 1 for different 7/t and t/I values. It can be seen that these models 
predict very different values of maximum stress. The shear stress difference can reach up to 40%, while the 
peel stress difference can be almost 30%. It appears that the differences between the two models are 
especially large for thicker adhesive layers and shorter overlaps. Interestingly, the sign change between t/ 
values of .001 and .01 in both plots indicates that there may be a t/] value which causes the difference to 
hover around zero. 
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Figure 2-6. Relative difference in maximum adhesive centerline (z,;=0) shear and peel stress between 
Models 3 and 2. 
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Figure 2-7. Comparison of the maximum shear and peel stress predicted by Models 2 and 1 for different 
joint configurations, along with points indicating the cases for comparison with 2-D FEMs. 


2-D plane stress FEMs were constructed for the four cases, each case with parameters as described in 
Table 2-1 and the points labeled in Figure 2-7. The results of the comparison between the FEM models and 
the joint element models are found in Figure 2-8 up to Figure 2-16. Figure 2-8 compares the reaction force 
predicted by the three adhesive models, normalized and compared with that of the 2-D element dense mesh 
model. Since the joints were given a prescribed displacement, this plot of the reaction force also represents 
how well the stiffness was captured by the joint models. This corresponds with the shear stress comparison 
discussed earlier, since the reaction force for a linearly elastic joint can be equated to the integral of the 
shear stress over the are in the area of the adhesive in the xy-plane. Therefore, model 1 deviated from the 
other models and the Abaqus model the most for cases 1 and 3, where the adhesive was the thickest. An 
investigation of the adhesive centerline stresses along with 2-D contour plots showing the shear and normal 
stress distribution in the adhesive as predicted by the FEM and three models for each case will be presented 
in the following section for each case for a more detailed discussion. 


2.3.1.1 Case 1 


Case 1 had a thick, long adhesive layer with 7/t = 1 and t/I = 0.001. The length of the adhesive displayed 
in Figure 2-9 and Figure 2-10 is only 0.5% of the length because the adhesive in the middle region is 
effectively stress free, which makes it this region unimportant for the current study. Obviously, this joint 
would be very inefficient because almost all of the stress is held by less than 1% of the adhesive. All three 
models predict the peel stress fairly well, although Models 3 and 2 appeared to be slightly more accurate. 
Looking at the adhesive peel stress distribution in Figure 2-10, Models 2 and 1 predict constant peel stress in 


the z-direction, while the peel stress distribution of the 2-D FEM varies nonlinearly in the z-direction, even 
past the free end. 
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Figure 2-8. Reaction force predicted by the three adhesive models compared to the 
2-D element dense mesh model for single lap joins with different geometries. 


e Abaqus 2-D elements 

--- Model 1 

4 + —-Model 2 ‘ 
—Model 3 eee A 


7; (MPa) 


0.495 0.496 0.497 0.498 0.499 0.500 0.495 0.496 0.497 0.498 0.499 0.500 
x/l x/l 


Figure 2-9. Adhesive centerline stresses predicted by different models for Case 1 parameters. 
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Figure 2-10. Adhesive stress contours predicted by different models for Case 1 parameters. 


Looking at the centerline shear stress, Model 1 does not seem to match up with the FEM very well. 
While Models 3 and 2 are low, they appear to capture the general trend. Obviously, none of the models 
reflect the free end condition of zero traction. The shear stress contour plots in Figure 2-10 show that Model 
1 has shear stress constant in the z-direction, but the FEM model shows that the stress varies significantly in 
the z-direction, which might explain the inaccuracy of the model. The distribution of shear stress in Models 
3 and 2 is a closer match to the FEM model, although the slopes of the lines of constant stress appear to be 
incorrect past the free end. 


2.3.1.2 Case 2 


Case 2 was composed of a joint with a thin, long adhesive layer (7/t = 0.1 and t/ = 0.001). The length of 
the adhesive displayed in Figure 2-11 and Figure 2-12 is only 0.5% of the length because only the ends 
display behavior worth comparison. All three models predict the peel and shear stress very well. The reason 
can be seen in the contour plots in Figure 2-12. Although it is difficult to see, the effect of the stress 
singularity and traction free boundary at the free end of the adhesive predicted by the 2-D FEM dies out 
very quickly, and the remainder of the adhesive distribution in the z-direction similar to that predicted by the 
joint element models. Since the thickness of the adhesive is so small, the stress can be effectively modeled 
as constant in the z-direction. Therefore, using the simplified Model 1 would still yield very accurate results 
for this case. 
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Figure 2-11. Adhesive centerline stresses predicted by different models for Case 2 parameters. 
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Figure 2-12. Adhesive stress contours predicted by different models for Case 2 parameters. 
2.3.1.3 Case 3 


Case 3 had a thick, short adhesive layer with 7/t = 1 and t/! = 0.1. The right half of the adhesive is 
shown in the plots. The three model predictions of the shear stress were not very consistent with the 2-D 
FEM. Model 1 vastly over predicts the stress levels, and Models 3 and 2 under predict the shear stress in 
most of the adhesive (Figure 2-13). Looking at the shear stress distribution of the FEM in Figure 2-14, the 
effects of the free end boundary do not die out quickly as with the thin adhesive cases. The constant z- 
direction stress distribution of Model 1 is insufficient to correctly model the adhesive, especially when it 
comes to the shear stress. 

The peel stress, on the other hand, was predicted reasonably well. All three models predicted similar 
trends for the adhesive centerline peel stress, even though the prediction appears to lag behind the stress 
predicted by the 2-D FEM. Generally, Model 1 would be a poor choice for a joint of this type. Models 2 and 
3 are closer, but still do not accurately predict the stress levels of the single lap joint for Case 3. 
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Figure 2-13. Adhesive centerline stresses predicted by different models for Case 3 parameters. 
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Figure 2-14. Adhesive stress contours predicted by different models for Case 3 parameters. 
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2.3.1.4 Case 4 


Case 4 was composed of a joint with a thin, short adhesive layer (7/t = 0.1 and t/] = 0.1). All three 
models under predicted the shear stress slightly, and surprisingly, Model 1 is the closest to the 2-D FEM 
(Figure 2-15). The contour plot of the FEM shear stress in Figure 2-16 shows that the effects of the free 
edge disappear quickly in the adhesive. The lines of constant shear stress appear linear, but Models 3 and 2 
have the wrong slopes, which probably accounts for the inaccuracy. Since Model 1 results in stresses 
constant in the z-direction, it is slightly more accurate than Models 3 and 2. The models predicted the FEM 
peel stress slightly more accurately, but there was still a lag similar to Cases 1 and 3. It is unclear what 
causes the lag; whether it is an effect created by the difference in free edge conditions, or simply due to 
inaccuracies in the model. 

The models presented here reflect the behavior of the joint in a global sense to be used for early sizing 
studies. When one keeps that in mind, the stress states predicted by the models for all of the cases appear to 
be adequate. The models are especially suited for thin adhesive layers, and appear to be even more accurate 
for long, thin adhesive layers. This tends to minimize edge and through-thickness effects, which are only 
approximated by the analytical models. 
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Figure 2-15. Adhesive centerline stresses predicted by different models for Case 4 parameters. 
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Figure 2-16. Adhesive stress contours predicted by different models for Case 4 parameters. 


NASA/CR—2012-217606 21 


2.3.2 Composite Adherends 


The composite adherend modeling capability of the joint element model was compared with the 
commercial sizing code HyperSizer [22], which uses an analytical formulation derived by Mortensen and 
Thompsen [6]. The two were compared for a bonded doubler under tension. The material properties are 
outlined in Table 2-2, the layup and boundary conditions can be found in Table 2-3, and the geometric 
parameters are shown in Figure 2-17. The joint element adherend and adhesive models chosen for the study 
were the Euler Bernoulli adherend model and adhesive model 1 (Winkler foundation), which are most 
consistent with the HyperSizer joint analysis method. Figure 2-18 shows a comparison of the normal 
stresses in the x-direction (0,) through the thickness of the upper adherend at a distance of half of the length 
of the stiffener from the supported edge. As can be seen, the stresses in the layers are nearly identical. 


Table 2-2. Material properties used in Hypersizer analysis [22]. 


E, E> E3 Gi G31 G23 
(GPa) (GPa) (GPa) (GPa) (GPa) (GPa) V12 V13 V23 


Boron/Epoxy 223 24.1 24.1 8.48 8.48 8.48 0.23. 0.25 0.32 
Epoxy 3.07 3.07 3.07 1.14 1.14 1.14 0.348 0.348 0.348 


Table 2-3. Layup and boundary conditions for Hypersizer analysis [22]. 


Loads and boundary 


Material Layup conditions 
1l 
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Figure 2-17. Geometric parameters for composite stiffened panel comparison with HyperSizer [22]. 
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Figure 2-18. Comparison of normal stresses in the x-direction of the mid-joint cross section 
of adherend 1 of a stiffened panel found using HyperSizer [22] and the joint element. 


2.3.3 Tapers and Steps 


In order to check the validity of the model and its ability to accurately predict stresses in a joint with 
tapers and steps, two joints were modeled with the joint elements and a dense mesh model using the 
commercial software Abaqus with 2-D plane stress elements. The joint element adherend and adhesive 
models chosen for the study was the Euler Bernoulli adherend model and adhesive model 1 (Winkler 
foundation).The first joint modeled was a stepped double-strap joint, which was chosen specifically to 
demonstrate how well the joint element could predict stresses in the vicinity of the steps. Second, a tapered 
joint was modeled to explore how well a smooth taper could be modeled using discrete, stepped joint 
elements. Stresses obtained using both models were compared for both joints, and conclusions were drawn 
in order to guide joint designers in the correct use of the joint element. 


2.3.3.1 Stepped Double-Strap Joint 


A stepped double-strap joint was modeled using reduced-order joint elements and a dense mesh with 
2-D plane strain elements. The geometric and material parameters used are found in Figure 2-19. The joint 
was modeled using 26 joint elements and around 99,000 2-D plane strain elements in Abaqus (Figure 2-20). 
The left end was clamped, and the right end was restrained from rotation and transverse displacement. The 
stresses in the upper and lower adhesive layers are shown Figure 2-21a and b, while the stresses in the upper 
and lower surfaces of the left and right adherends are shown in Figure 2-21c and d. 
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Figure 2-19. Geometric and material parameters of the stepped double-strap joint. 
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Figure 2-20. Two models compared in the study, (a) joint element model using 26 elements and (b) Abaqus 2-D 
plane stress model using 99,000 elements. 


The shear and normal stresses in the centerline of the adhesive layers agreed very well between the two 
models, even in the regions near the steps. The only discrepancies were found near the edges of the 
adhesive. The 2-D Abaqus model treats the adhesive as a continuum, while the joint element treats the 
adhesive as normal and shear springs. Therefore, stress singularities and traction free boundary conditions 
are not captured in the joint element model. However, for long overlap regions and thin adhesive layers, the 
difference is minimal. 

The adherend stresses are also in good agreement (Figure 2-21c and d) except at the reentrant corners 
and the loaded ends. As stated above, the stress singularity of the reentrant corners is not captured by the 
model. However, overall trends are captured (reaction force for the Abaqus 2-D model was 4.323 KN, while 
the joint element model predicted 4.326 kN), making this model a good tool for initial design and sizing, 
particularly considering that in reality, the spew fillet would mute the discrepancy at the reentrant comers. 

It is the stress in the doublers which requires special attention. First, consider the surface of the doublers 
away from the surface of the steps (Figure 2-22a and d). The stress peaks in the joint model around the steps 
appear to correlate quite well with the stress predicted by the Abaqus 2-D model. Near the steps, the joint 
element model, which treats the adherends as beams, predicts stress oscillations. The FE model, which treats 
the adherends as continua, does not predict these oscillations at this surface, which is some distance from the 
steps. 

The stress on the surface of the steps (Figure 2-22b and c) shows fairly good correlation between the 
models, except for the last step (the outer peaks in Figure 2-22b and C). It is believed that since the joint 
element model does not take the traction free surface of the step into account, the difference between the 
two models grows as a larger percentage of the thickness is a stress free surface. At the last step, half of the 
surface is traction free, which causes a discrepancy between the two models. However, despite the fact that 
the joint element model doesn’t capture stress singularities at reentrant comers and ignores traction free 
surfaces of the steps, the joint element model appears to correlate well with dense 2-D finite element 
meshes. 
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Figure 2-21. Stresses in the (a) upper and (b) lower adhesive centerlines, and (c) left and 
(d) right adherends in the double-strap joint. 
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Figure 2-22. Stresses in the (a) lower and (b) upper surfaces of the top doubler, and (c) lower and (d) upper surface 
of the bottom doubler in the double-strap joint. 
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2.3.3.2 Tapered Joint 


A tapered joint was also modeled with both models to see how well the smooth taper could be 
approximated by a discrete stepped model using joint elements. The parameters for the joint are illustrated 
in Figure 2-23. The joint was modeled using 7, 12, and 22 joint elements (Figure 2-24a-c). This was to show 
the effect of mesh density in the vicinity of the joint. All three were compared with a 120,000 2-D plane 
stress element model in Abaqus (Figure 2-24d). 

The stresses in the upper and lower surface of the left adherend near the taper are shown in 
Figure 2-25a and b, while those of the right adherend can be found in Figure 2-25c and d. The normal and 
shear stresses in the adhesive layer are shown in Figure 2-25e and f, respectively. The stresses in the stepped 
joint model appear to oscillate around those of the smooth taper Abaqus 2-D model in the adherends 
(Figure 2-25a-d). However, as the density of the joint elements increases, the stress predicted appears to 
converge on the stress in the Abaqus 2-D smooth taper model. 
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Figure 2-23. Geometric and material parameters of the tapered joint. 


Joint Element 
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(b) 12 Elements 
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Figure 2-24. Models compared in the study: joint element using (a) 7, (b) 12, (c) and 22 elements, and 
(d) Abaqus 2-D plane stress model using 120,000 elements. 
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Figure 2-25. Stresses in the (a) left and (b) right adherends and (c) peel and 
(d) shear stress in the adhesive of the tapered joint. 
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The peel stresses in the adhesive centerline appear to be quite different between the Abaqus 2-D model 
with a smooth taper and the stepped joint element models. Both have maximum stresses at the ends of the 
adhesive layer, but the joint element model shows a near zero stress away from the ends while the Abaqus 
2-D model has a nearly constant stress of about 8 MPa in the middle of the joint. This is due to the way the 
joint element represents the smooth taper. The adhesive layer is modeled as horizontal discrete pieces of 
adhesive only, while the vertical component is ignored. The vertical component would carry a near-constant 
peel stress, which is why the peel stress in the middle of the joint is not captured. Therefore, one must keep 
in mind that this component of stress will not be captured by the joint element model. However, the joint 
element model appears to converge to the Abaqus 2-D solution for the shear stress in the adhesive layer. 


2.4 Conclusions 


Eccentric load paths caused by single lap joints often require a detailed, dense mesh in structural finite 
element models of assemblies containing single-lap joints, which can be costly for global, vehicle scale 
models. To reduce computational time required for these large-scale models and aid joint sizing early on in 
the design phase, a joint element was created to model the behavior of a joint with a single element. This 
was accomplished by embedding an approximate analytical solution of the single lap joint into a finite 
element. The model discussed in this chapter assumes that all materials remain linearly elastic, but the 
extension of this model to include progressive failure discussed in Chapter 4 will require an iterative 
solution procedure. With this in mind, it is imperative that the simplest formulation be used to reduce the 
complexity and computing time required for each iteration. Therefore, three adhesive models of varying 
complexity were compared. 

The first adhesive model, Model 1, treats the adhesive as a bed of uncoupled linear springs (Winkler 
foundations), making the stresses independent of the z-direction. Model 2 couples the normal springs 
together so that the variation of w displacements across the joint aids in shear strain accumulation. Finally, 
Model 3 assumes that the adhesive is in plane stress, and considers gx, Gazz, ANd Taxz. Additionally, a model 
assuming plane strain was created, but not compared in this study. All models assumed the displacements in 
the adhesive layer to vary linearly in the z-direction. The maximum peel and shear stress predicted by these 
three models was compared to ascertain how much of a difference each additional assumption matters. It 
was found that there was little difference between Models 2 and 3. Since Model 2 is simpler than Model 3, it 
can be concluded that Model 3 is unnecessary, unless the value of o,,, is desired. On the other hand, Models 
1 and 2 predicted very different maximum stresses, especially for joints with thick adhesive layers and short 
joint overlap lengths. 

Four joints with different combinations of thick and thin adhesive layers and long and short overlap 
lengths were studied in depth and compared with 2-D finite element models to 1) identify which models 
were more accurate, and under what parameter conditions and 2) to show that the joint finite element could 
be used to give a general prediction of joint behavior. It was found that all of the models were more accurate 
for joints with thin adhesive layers. Also, all of the models were slightly more accurate for longer joint 
overlaps. Therefore, for thin joints, Model 1 should be used because it can predict an answer just as accurate 
as the other models, with a much simpler formulation. For thick joints, Model 2 would be preferred because 
it is more accurate than Model 1, although the models are generally not as accurate for thick joints. 
However, most joints in application have very thin adhesive layers, so Model 1 should be adequate for most 
real-life situations. It was also shown that for all cases, the joint elements were more than adequate at 
predicting the behavior of a joint for early design sizing purposes. 

A second validating example compared the joint element with HyperSizer [22] sizing software for a 
stiffened composite panel. The comparison showed an almost perfect match of the stresses in the upper 
adherend layer. 
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Finally, it was shown how well the joint element can model joints with tapers or steps. For discrete 
steps, the joint element appeared to agree with the 2-D dense mesh finite element solution. For a gradual 
taper, it was illustrated that the stresses in the adherends predicted by the joint element model appear to 
oscillate around the dense mesh solution, and that the stress converges to the solution by increasing the 
number of steps. Also, the lower bound of the adhesive shear stress is a good representation of the smooth 
taper solution. The average peel stress is not captured using the joint element model, but the maximum 
stress at the ends of the adhesive are still predicted, and these high stresses are of concern in predictions 
related to joint strength and durability. 
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Chapter 3 
Functionally Graded Adhesives for Composite Joints 


Adhesives with functionally graded material properties are being considered for use in adhesively 
bonded joints to reduce the peel stress concentrations located near adherend discontinuities. Several 
practical concerns impede the actual use of such adhesives. These include increased manufacturing 
complications, alterations to the intended grading due to adhesive flow during manufacturing, and whether 
changing the loading conditions significantly impact the effectiveness of the grading. An analytical study is 
conducted to address these three concerns. An enhanced joint finite element, which uses an analytical 
formulation to obtain shape functions, is used to model the joint. Furthermore, proof—of-concept testing is 
conducted to show the potential advantages of functionally graded adhesives. In this study, grading is 
achieved by strategically placing glass beads within the adhesive layer at different densities along the joint. 


3.1 Introduction 


One major drawback of adhesively bonded joints is that the load path eccentricity causes the appearance 
of peel stress concentrations at the end of the adhesive layer. There has been a vast amount of research 
conducted in an attempt to reduce these stress concentrations, such as tapering the end of the adherend [23], 
increasing thickness of the adhesive at the end [24], fillets [20], novel joint geometries [25], and joint 
insertions [26], to name a few. All of these methods involve local details of adherend geometry (except the 
adhesive fillets), which typically increases part complexity. 

Material grading occurs in nature at material interfaces to reduce stress concentrations [27]. Biological 
interfaces such as tendon to bone joints have been found to have graded material properties to distribute 
stress more evenly across the joint [28]. In this same spirit, material grading has been applied to 
adhesively bonded joints. Although grading the adherends has shown promise [29], many more researchers 
have investigated grading the adhesive properties. Some of the earliest grading of the adhesive was reported 
by Patrick [30] and Raphael [31], where grading was discretely achieved using two adhesive materials (i.e., 
bi-adhesive). Recently many other researchers have investigated bi-adhesive joints with mixed results. 
Sancaktar and Kumar [32] graded the adhesive by making rubber toughened regions, and found that 
the selectively toughened joints had the same strength as the fully toughened joints. Piresa et al. [33] used 
two adhesives to bond aluminum, and found up to a 22% increase in joint strength. Fitton and Broughton 
[34] bonded carbon reinforced plastic (CFRP) to steel, and found that it was crucial to optimize the amount 
of each adhesive and that some configurations did not benefit from grading. Da Silva and Lopes [35] looked 
at the influence of the ductility of adhesives on bi-adhesive joint strength experimentally and theoretically. It 
was found that the bi-adhesive joints out-performed joints with the more brittle adhesive alone, but did not 
always improve on the more compliant single adhesive joints. More recently, Kumar and Pandey [36] 
modeled bi-adhesive joints using nonlinear 3D finite elements compared with a 2D finite element model, 
and found that the 2-D model could not fully capture the complex multi-axial stress state. Valleé et al. [37] 
investigated bi-adhesive joints, among other stress reduction methods, and found that the adhesive stress 
was not linked to the joint strength of the configurations tested (which displayed adherend failure). 

It appears that the first instance of grading the adhesive with a non-stepwise function was Kumar [38]. 
This purely theoretical investigation first compared a continuous (non step-wise) functionally graded 
adhesive (FGA), where the modulus was graded using a quadratic function, with a step-wise graded 
equivalent for different overlap lengths and adhesive thicknesses in a tubular joint. It was found that the 
continuous FGA reduced the shear and peel stresses in all cases. Second, four “arbitrarily chosen” functions 
were compared to show that the grading function can be manipulated to optimize joint performance. The 
current study aims to increase the understanding of both bi-adhesive and non-stepwise FGAs to make them 
a more viable, realistic, and advantageous choice for actual application in composite structures. 
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Some potential drawbacks to FGAs have also been discussed in the literature. The first two were 
identified by Hart-Smith [24], where he pointed out that: 1) small gains over just using the ductile adhesive 
alone may be inadequate when considering the production difficulties and 2) during manufacturing, there is 
the “inevitable tendency for the stiff adhesive to squeeze out and displace the ductile adhesive,” making it 
probable that the resulting joint will be worse-off than using the ductile adhesive alone. A third concern was 
raised by Aboudi et al. [39] while investigating the response of metal matrix composites with tailored 
microstructures. They found that functionally grading the properties of a material may be detrimental when 
the loading is changed, such that the stress gradient in the material is reversed; this is a valid concern for 
practical situations where all loading cases are not necessarily known. 

The current study aims to address these practical concerns to show that FGAs are a viable means of 
decreasing the peel stress in an adhesively bonded joint. An analytical model is constructed and used to 
compare the stresses in a butt-end joint configuration with four different functions of modulus graded 
adhesive: constant (single adhesive), discrete (bi-adhesive), linear, and exponential. First, the potential gains 
in stress reduction for FGA joints over joints with a single adhesive are shown. Along these lines, it is 
shown that additional stress reductions can be achieved by lowering the modulus of the more compliant 
adhesive. Since it is likely that step-wise grading will appeal from a manufacturing perspective, stress 
reduction of a step-wise graded adhesive with many steps is investigated (with single or bi-adhesive being a 
special case). Second, the study addresses the issue of adhesive flow during bonding by showing how 
sensitive the optimum for the three FGAs is to perturbation of the grading. Third, multiple load scenarios are 
examined to address the concern of changing loading conditions. Results indicate that the stress magnitude 
gradient in the adhesive remains unchanged under different loading conditions, making joints a perfect 
application for material grading. Additionally, it is shown that the grading can still be optimal under 
different loading cases. Addressing these three concerns provides significant impetus for the use of FGAs in 
industrial applications. 

The model used to obtain the adhesive stress for different FGAs is a structural finite element made 
specifically for adhesively bonded joints. Motivated by the desire to create a computationally efficient tool 
for designing joints within a coarse, vehicle scale finite element model [13], we combined an analytical 
formulation with a finite element in the joint element. This concept has been often referred to as the exact 
stiffness matrix method, and has been previously applied to the beam on an elastic foundation problem 
[11,12]. The joint element is capable of capturing the stresses in a mesh-independent, efficient manner. Such 
an efficient method is pivotal to the current study, allowing the analysis of over 20,000 different joints for 
optimization and parametric studies on a desktop computer in a fraction of a second per joint. The 
formulation adopted here is altered from the previous section to account for a graded adhesive modulus and 
is presented below. A linear elastic material model is used for several reasons, simplicity being the most 
prominent. Also, since a controlled method and material system for manufacturing FGAs has not yet been 
identified, failure of the joint and post-failure response is not addressed. Finally, it should be noted that after 
initial departure from material linearity (due to damage or plasticity) and before crack formation, the 
adhesive modulus is effectively a continuous function across the joint, which causes more load to be 
transferred to the inner regions of the joint. However, the main idea of a FGA is that this effect can be 
achieved without taking on irrecoverable damage. Since the benefit of FGAs can be realized without 
material damage, this study will be limited to the linearly elastic regime of the adhesive. Geometric 
nonlinearity is also ignored for simplicity and because it is not expected to have a large effect on a stress 
comparison study between adhesive systems. 

The analytical findings are complemented with an experimental “proof-of-concept” testing to illustrate 
the benefits of FGAs. The adhesive was graded by adding different volume percentages of glass beads, 
although no precise method was used to control the grading except the eye and hand of the person preparing 
the specimens. By showing that a joint can benefit from grading in such a rudimentary manner, the potential 
for more drastic gains through controlled and precise grading can be argued for. 
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3.2 Method 


3.2.1 Formulation 


The formulation for the joint element with a functionally graded adhesive is almost identical to the 
formulation of the element with a non-graded adhesive as shown in Section 2.2. Therefore, only details that 
diverge from the previous derivation will be shown here. 

The adhesive and adherend models utilized for this study are the Euler-Bernoulli adherend model 
(Section 2.2.1.1) and the simple Winkler Foundation adhesive model (Section 2.2.2.1), although any 
combination of adhesive and adherend models could be used. The adhesive moduli are now considered to be 
functions of x, causing the stresses in the adhesive to be 


O,, = Ey, (Ez, 3-1 


and 


T= GAY, 3-2 


This causes the system of governing equations to now have non-constant coefficients (compare with 
Equation 2-32) 


u, =A(x)u oe 


where the coefficient matrix, A(x), is now a function of x. Since the coefficient matrix is non-constant, a 
new method is needed to solve the system of ordinary differential equations. 


3.2.1.1 Method of constant segments for solving linear, homogeneous, non-constant coefficient system 
of ordinary differential equations 


In order to solve the system of ordinary differential equations, a semi-numerical method of solution was 
adopted. Traditional differential equation solving techniques employing numerical boundary conditions 
could not be employed because the boundaries (nodes) contain unknown, symbolic conditions. Therefore, 
the domain was split into segments in which the coefficient matrix, A(x), is considered constant and solved 
using the matrix exponential. First, consider segment n+1, with a local x-direction coordinate system x’ 
which originates at the left side of the segment, x=x, (Figure 3-1). The other end of the segment is at x ’=Ax 
and x=Xp+1. 


Figure 3-1. Method of linear segments involves breaking up a function into small segments 
and assuming the function is constant within the segment. 
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It is assumed that Ax is sufficiently small so that A(x) can be considered constant within each such 
segment. The linearized coefficient matrix A,.;, is taken to be the coefficient matrix evaluated at the 
midpoint of the segment: 


Avi = A(x, as =e) . a 


Within segment n+1, the system can now be expressed as a system of ordinary constant coefficient 
differential equations in the local coordinate system, x’, of the form: 


3-5 


As before, there are 6N eigenvalues of A,.1: N real eigenvalues, 2N complex eigenvalues, and 3N 
repeating eigenvalues. Therefore, the solution is made up of N exponential terms, 2N exponential terms 
multiplied by a sine or cosine, and the 3N repeating eigenvalues correspond to a third order polynomial 
found in a standard beam solution. Such a complex solution shows that merely employing standard beam 
shape functions to the joint problem would be inadequate in capturing the nature of the complete solution. 


The solution of the system in Equation 3-3 can be written in terms of the matrix exponential, e*""”’ , and 
a vector of unknown constants, C,+1, as 


u(x')= ia Seer i 
The matrix exponential can be expressed as the infinite series [17] 
oo xk 3-7 
Anax = y—a‘ 
x20 K! 


To get faster convergence, a method of scaling and squaring [18] is employed, and the series is calculated 
up to a value of k which yields an acceptable error, € . The error can be defined many ways, but the current 


study defined the error as the difference between the 1-norms of e*"’ for k-1 and k. The value of the 
acceptable error was set at € = 0.0001. 
In the local coordinate system, the solution u, at x’=0 can be expressed as 


3-8 
u, = ea Oper = Gaui ; 
and the solution, u,+:, at the end of the segment (x ’=4x) can be written in terms of the solution at the 
beginning of the segment, 
AnytAx 3-9 


to eliminate the vector of constants, C,,.;. Similarly, the solution at the end of the previous segment can be 
expressed as 


- 3-10 


and so on, down to the first segment, which has the solution: 
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u, =e"C, 3-11 
Therefore, the solution at any segment, n+1, can be expressed in terms of the vector of constants from the 
first segment, Co, by the equation: 

= e OSG: 3-12 


Una 


where 


n+1 3-13 
eNoinet TT ean 
m=0 


To demonstate this method, Figure 3-2 shows a convergence study for the differential equation 
y= (x° —e* sin 4x) y with initial conditions y(0) =1 and y,,(0) = 0. It can be seen that the percent 


error between the method of constant segments and the solution found using a 4" order Runge-Kutta [40] 
solver with a step size of 0.001 diminishes with the number of segments. For the current study, 200 
segments were used in each element. 

It should be acknowledged that the discretization involved in this method is very similar to conventional 
finite element discretization. This may serve to reduce the advantages of using the joint element, but it still 
allows for the joint to be represented with one element, simplifying analysis steps like mesh generation and 
post-processing. 
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Figure 3-2. A comparison of solutions to the homogeneous, linear, non-constant coefficient differential equation 
Vx= (x3 —e* sin 4x) y with initial conditions y(0)=1 and y ,(0)=0 (a) using the method of constant segments 


with two, four, and six segments and (b) the relative error as a function of the number of segments used. 
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3.3 Results and Discussion 


3.3.1 Baseline Configuration 


To assess the performance of functionally graded adhesives (FGAs), a baseline joint configuration was 
identified based on the configuration in the proof-of-concept testing (Section 3.3.5). A single strap joint (or 
butt end joint) was chosen because of its ease of manufacturing, symmetry, and single dominant stress 
concentration in the middle of the joint. The geometric and material properties are defined in Figure 3-3a, 
and the finite element representation is shown in Figure 3-3b. Only half of the joint was modeled due to 
symmetry, and the overlap section was modeled by one or several joint elements (depending on the 

number of discrete regions of continuous adhesive modulus) while the non-joint adherend section is 
modeled with one beam element. The loaded end is constrained from rotation and vertical translation, while 
the symmetric face of the doubler is constrained from horizontal translation and rotation. 

The values of the material and geometric parameters used for the analytical modeling are found in 
Table 3-1, and are based on realistic values that coincide with the proof-of-concept test specimens discussed 
in the next section. The FGAs were compared with two different single adhesive systems (Figure 3-4a). 
These two adhesives provided upper and lower bounds for the grading functions, and will be referred to as 
E,, and E,, respectively. The grading functions chosen for investigation included a step (bi-adhesive), linear, 
and exponential function. These functions were reduced to a single grading variable, |, as defined in 
Figure 3-4b, c, and d. For the step function, | is the length of the more compliant adhesive. For the linear 
function adhesive, / is the length of the section in which the modulus decays linearly, and for the exponential 
function it is the length of adhesive which has a modulus less than 99% of E,,. Although these may not be 
the optimal grading functions, these functions were chosen because of their ability to be reduced to one 
variable, allowing for simple and clear sensitivity studies. It was assumed that regardless the grading, the 
relationship between the Young’s modulus and shear modulus remained constant or, in other words, that the 
Poisson’s ratio remained constant. A similar assumption was made by Apetre et al. [41] for functionally 
graded sandwich beam cores. 
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Figure 3-3. Single strap joint (a) geometric and material parameters and 
(b) joint finite element representation assuming symmetry. 


Table 3-1. Parameters of the baseline single strap joint configuration used for the theoretical study. 


P la lo t b n E Ey E, v Tnet 
(kN) (mm) (mm) (mm) (mm) (mm) (GPa) (GPa) (GPa) (MPa) 
4000 82.6 38.1 1.1 25.4 0.4 108.5 2.5 i 0.34 8.3 
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Figure 3-4. Single strap joints with (a) constant modulus adhesives were compared with joints with functionally 
graded adhesives, including (b) step-wise graded, (c) linearly graded, and (d) exponentially graded. 


Most stress values reported in the theoretical study were normalized by the net shear stress, Tne, defined 
by: 


P 
Tet = Lb 3-14 


3.3.2 Stress Comparison 


The first goal of this study was to address the concern that functionally graded adhesives (FGAs) are not 
worth the increase in manufacturing difficulty. This section addresses this concern by showing the potential 
adhesive peel stress reduction due to FGAs over single adhesive joints. This will be shown for the three 
functions defined previously (Figure 3-4). Next, the result of changing the lower bound modulus, E;, will 
show how further stress reductions can be achieved. Finally, since stepped adhesives will most likely be the 
easiest FGAs to manufacture, the effect of increasing the number of steps and spacing optimization will be 
investigated. All optimization was carried out using a genetic algorithm. 

The grading parameter, |, was optimized to reduce the maximum peel stress in the adhesive for the three 
FGAs, and the resulting moduli are plotted in Figure 3-5a. As shown, the region of gradation is very small, 
about 2% of the overlap length for the step and linear function adhesives, and around 5% for the 
exponential. The peel stress in the adhesive for half of the symmetric joint is plotted in Figure 3-5b for each 
adhesive. The single adhesive joints are in blue, and the FGAs are in black. The step FGA has two stress 
peaks; one at the end of the adhesive and one at the interface between the two adhesives. The linear and 
exponential FGAs have a rounded stress peak, and appear to result in very similar stress distributions. 
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Figure 3-5. Optimized configurations for the single strap joints for different functions of graded adhesive compared: 
(a) modulus across the adhesive and (b) centerline peel stress across the adhesive. 


The maximum stress in the adherend, doubler, and adhesive for the single adhesive joints and the 
optimized FGA joints is found in Figure 3-6. All stress values are normalized by the maximum stress found 
in the stiffer single adhesive for viewing all stress components in the same plot. The stress reported for the 
adherend and the doubler is the normal stress in the x-direction, and the maximum value of the stress is 
found at the upper and lower surfaces of the adherend and doubler, respectively. With composite laminate 
adherends, the most important stress in the adherend is usually the peel stress (z-direction) between the plies 
because this is where failure often initiates. However, the current configuration contains only one ply, so 
there is no such interlaminar shear stress. 

Some important aspects of using FGAs are illustrated in Figure 3-6. First, the FGAs in this study 
outperformed the stiffer single adhesive joint, E;. Adhesive stresses were considerably lower and adherend 
and doubler stresses were not significantly impacted. This is important because the more compliant single 
adhesive joint had lower adhesive stresses, but higher adherend stress. The FGAs were able to lower the 
adhesive stress without affecting the adherend stress. Second, when compared with the E, adhesive, the 
FGAs reduced the adhesive peel stress but did not necessarily reduce the shear stress. It will be shown later 
that the optimum value of the grading variable | is not the same for minimizing peel as it is for minimizing 
shear stress. Therefore, the relative levels of peel and shear must be considered when designing FGAs so 
that the dominant stress can be minimized. However, typical adhesives are more ductile under shear loading, 
so peel stresses will normally be the minimized variable. 
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Figure 3-6. Comparison of the maximum stresses in the joints with different adhesives, where all maximum stress 
values are normalized by the maximum stress value found in the stiffer constant modulus (single adhesive) joint, 
and adherend and doubler stresses refer to those in the x-direction. 
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Figure 3-7. The maximum stresses in joints with linear FGAs and different lower bound modulus, El, where all 
maximum stress values are normalized by the maximum stress value found in the stiffer constant modulus (single 
adhesive) joint, and adherend and doubler stresses refer to those in the x-direction. 


In order to further reduce the maximum adhesive peel stress, the lower bound modulus can be 
decreased. The affect of changing the lower bound modulus on the maximum joint stresses is shown in 
Figure 3-7 for the linear FGA only (Figure 3-4c). For each value of E;, the grading parameter | has been 
selected to produce for the smallest maximum peel stress. By grading the adhesive from a very low modulus 
to a higher one, the peel stress can be lowered significantly, with peel stress reductions in this case of up to 
90%. At the same time, the shear stress is also significantly lowered. For low lower bound modulus values, 
the maximum adherend and doubler stresses are also impacted. The maximum adherend stress increases 
with decreasing FE), while the maximum doubler stress decreases slightly. While it may appear that the 
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Figure 3-8. Increasing the number of adhesives used for a discretely graded adhesive and optimizing the 
adhesive spacing lowers the maximum adhesive peel stress; (a) stress vs. number of adhesives, and 
(b) modulus at the end of the adhesive for 20 steps compared with the other functions. 


adherend stress is raised considerably while the doubler is only slightly decreased, it should be kept in mind 
that these values are normalized by the stresses found in the stiffer single adhesive system. The maximum 
doubler stress is around three times higher than the maximum adherend stress, so the stress is being taken 
from the highly loaded doubler to the less loaded adherend, which is advantageous when they are all made 
of the same material as in this case. 

Since it is likely that grading of adhesives will be accomplished using a step function with different 
adhesives for each step, the effect of the number of steps used on the adhesive peel stress is plotted in 
Figure 3-8a. For each point, it was assumed that the moduli of the steps are equispaced between E; and Ey. 
The length of each step was optimized using a genetic optimization routine to reduce the peel stress. Using 
more steps led to a lower maximum peel stress, eventually even lower than the linear and exponential 
functions. The optimized moduli for 20 steps is shown in Figure 3-8b, along with the optimal linear and 
exponential moduli. This sheds a little light on a better grading function, consisting of a constant modulus 
region followed by an exponentially higher order declining region. 

This section has shown that there are significant stress reductions possible by using FGAs. Unlike 
simply using a lower modulus adhesive, FGAs do not have a significant effect on adherend and doubler 
stresses. Furthermore, lowering the lower bound modulus leads to even greater stress reductions. These 
stress reductions are possible using stepped modulus adhesives, although the length of each step needs to be 
optimized. 


3.3.3 Sensitivity Study 


The second goal of this study was to address the concern that during manufacturing, the adhesive is 
pressurized and heated, often causing the adhesive to flow and even squeeze out of the joint. If a 
functionally graded adhesive (FGA) is specifically designed for a certain joint, such squeeze-out could 
either change the shape of the grading, the lower bound modulus, or both. This could result in an FGA 
which has higher stress than using the more compliant adhesive alone. This section seeks to address this 
concern by presenting a grading sensitivity study. For the purposes of this study, sensitivity will refer to the 
narrowness of the range of | values that results in lower maximum stresses than those obtained using the 
lower bound adhesive, E,. First, the effects of changing the grading parameter / on the maximum adhesive 
stress are quantified (Figure 3-9), followed by an investigation into changing the lower bound modulus E, 
and | (Figure 3-10). Through this study, insight is gained into the sensitivity of the grading and which 
functions are more tolerant to changes in the grading. 
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Figure 3-9. Observing the maximum adhesive (a) peel and (b) shear stress as a function of | as defined in Figure 3-4 
shows how sensitive the maximum stress is to the shape of the grading. 


The effects of changing the grading parameter / on the maximum adhesive peel and shear stresses for 
the three FGAs are found in Figure 3-9a and b respectively. For all three FGAs, the maximum shear stress is 
less sensitive to the grading parameter than the peel stress. Also, the optimum value of I is always greater for 
minimizing the shear stress than the peel stress. This is most likely because the peel stress peak is much 
more concentrated than the shear stress peak, so a steeper gradation is needed to minimize the peak. Also, if 
| becomes too short, the stress goes above the stress that would be found in a joint with just the more 
compliant adhesive. Unfortunately, with the addition of pressure, a decreasing / is more likely. The step 
function adhesive was the most sensitive: only a very small range of values of / results in lower stresses than 
just using the more compliant adhesive, E,. The linear function adhesive was not as sensitive, and it 
converges to a stress less than that of E; when | is large. Finally, the exponential function adhesive had a 
broad range of | values resulting in low maximum stress, making it the most tolerant to changing the grading 
parameter. 

When the joint is being manufactured, the adhesive is heated and pressure is applied. The most likely 
result will be that the adhesive will spew out and the modulus of the adhesive at where the adherend ends 
will be greater. At the same time, the grading parameter / may also change. Therefore, it is important to 
know how sensitive the stress is to increasing the modulus and how that affects the grading parameter 
sensitivity. Based on Figure 3-10, the sensitivity of the maximum adhesive stress for all three FGAs to | 
increases with increasing lower bound modulus E;. Additionally, the optimum grading parameter ! decreases 
as the lower bound modulus increases. Therefore, it is recommended that while using FGAs in joints, one 
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Figure 3-10. The sensitivity of the maximum adhesive peel and shear stress to decreasing the lower bound modulus, 
E,, for the discrete (a and b), linear (c and d), and exponential (e and f) FGAs. 


should either figure out how much flow will occur and plan accordingly, or use a grading system which will 
not move when adhesive flows. An example of such a system is the glass beads used for the proof-of- 
concept testing which will be described in Section 3.3.5. Since beads are used for thickness control also, 
applying pressure causes the beads to be clamped between the adherend and doubler and unable to move. 
Such a system of grading is nearly independent of adhesive flow. 


3.3.4 Effect of Loading 


The third goal of this study was to show the effect of changing the loading conditions when using FGAs 
in joints. In many applications of functionally graded materials, a change in loading can lessen or even 
reverse the stress gradient, rendering the material grading useless or even detrimental. To show that this is 
not the case with adhesively bonded joints, six different loading scenarios and their respective adhesive peel 
stress distributions are shown in Figure 3-11. The boundary conditions were altered from that of the baseline 
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Figure 3-11. For different loading scenarios at the end of the joint ( (a) axial load, (c) moment, and 
(e) vertical load) the stress gradients of the adhesive peel stress (b, d, and f) remains the same direction; 
high at the ends and low in the middle. 


configuration (Figure 3-3) to represent an end moment and shear force. The end of the doubler was 
clamped, with no other constraints. The peel stress in the adhesive layer for a tensile and compressive axial 
load, positive and negative moment, and up and down vertical shear force is shown in Figure 3-11b, d, and f 
respectively, while the loading, deformed configuration, and boundary conditions are shown in 

Figure 3-11a, c, and e. 

For all of these loading cases, the gradient of stress goes from high stress (tensile or compressive) at the 
ends to low stress in the middle. Therefore, functionally grading the adhesive in the manner prescribed will 
always be beneficial under mechanical loading; thermal loading conditions will be addressed in the future. 

Compressive peel stresses are not considered to initiate or propagate cracks in an adhesive layer. 
Therefore, when the loading results in a compressive peel stress at the end of the joint, the subsequent stress 
“trough” (see Figure 3-5 for an example of a stress distribution and trough), which is comprised of tensile 
stress, becomes the critical stress. Although this trough is always significantly lower in magnitude than the 
peak stress at the end of the joint, the gradation should be designed to minimize the trough. By minimizing 
the trough rather than the peak stress at the end of the joint, the strength of a joint with a compressive peak 
stress at the end of the adhesive can be maximized. 

To simplify the parametric studies, the FGAs investigated were only designed to minimize the stress 
concentration in the middle of the joint (x/l,=0.5) because this was the dominant stress concentration in the 
baseline configuration. However, Figure 3-11d and f shows that changing the loading can cause both ends to 
be highly stressed. Therefore, it is highly recommended that the grading be applied to every area of a joint 
where there will be a stress concentration; i.e. at the ends of the adherends and doubler. Applying grading to 
all ends of the joint will ensure that the FGA will still be optimal regardless of loading conditions. 

Although Figure 3-11 does not explore all possible loading conditions, it is believed that these represent 
a large and varied enough sample to draw general conclusions. Different loading conditions away from the 
actual joint overlap region can always be resolved into a shear load, axial load, and moment near the joint 
overlap as has been done for many classical analytical formulations [21]. Therefore, changing the loading at 
the end of the adherend results in merely changing the relative magnitudes of the load components near the 
joint overlap region. 
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Figure 3-12. Maximum peel and shear adhesive stress and its sensitivity to grading parameter, 1, for different 
amounts of end moment (a, b) and end shear force (c, d) for joints with a linear FGA. 


To study the effect of changing the loading on the optimal grading shape and grading sensitivity, the 
normalized maximum normal (Figure 3-12a, c) and shear (Figure 3-12b, d) stress in the adhesive for joints 
with the linear FGA with different grading shapes was plotted. Figure 3-12a and b show the effect of 
different ratios of an end vertical shear force to axial force, and Figure 3-12c and d show the effect of 
different ratios of end moment to axial force. The loading was chosen so that the maximum stress in the 
adhesive remained at x/l,=0.5. As can be seen, the loading investigated here did not affect either the 
optimum grading or the sensitivity of the grading shape. Therefore, it can be concluded that the same 
grading would be beneficial for multiple loading scenarios. 

In addition to varying the loading scheme, the joint type was also altered to show that the adhesive peel 
stress magnitude gradient always goes from high stress near the adherend or doubler discontinuities to low 
stress elsewhere. The normalized adhesive peel stress for a single lap joint, bonded doubler, and double 
strap joint for different loadings is shown in Figure 3-13a, b, and c, respectively. This figure not only shows 
that the adherend peel stress magnitude gradient remains in the same direction for different joint types, but 
also suggests that the normalized stress gradient is very similar to that of the single strap joint. It might be 
the case that the same adhesive gradient could be utilized for not only different loading conditions, but 
different joint types. If it was the case that one grading type was beneficial for many different joint types and 
geometries, graded adhesives could be mass produced for general application rather than custom designed 
for each individual joint type and geometry. 


3.3.5 Proof-of-Concept Testing 


To illustrate the potential gains of FGAs, proof-of-concept testing was performed on single strap joints 
(the baseline case). The adherends and doubler consisted of a single ply of 0/+45 triaxially braided 
composites, with the axial tows consisting of around 80,000 fibers, biased tows with 12,000 fibers, a fiber 
volume fraction of 52%, and matrix of Epon 862 epoxy resin [42]. The axial 0° direction was aligned with 
the global x-direction as defined in Figure 3-3. The adhesive used was AF 163-2k [43] and the bond line 
was two layers thick in order for glass beads, normally used for thickness control, to be inserted in between 
(Figure 3-14). 
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Figure 3-13. For different loading scenarios for a (a) single lap joint, (b) bonded doubler, and 
(c) double strap joint, the stress gradients of the adhesive peel stress (b, d, and f) 
remains the same direction; high at the ends and low in the middle. 
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Figure 3-14. Diagram of single strap joint with placed glass beads to test functionally graded adhesive concept. 
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All specimens were manufactured together as a plate, with individual specimens cut out of the plate by 
means of a wet saw with a diamond-coated blade. The adherends and doubler were pre-infiltrated and cured 
prior to bonding. After bonding, the adhesive was cured in a heated press at 177° C under 0.43 MPa for 120 
minutes, and specimens were not handled for 24 hours. The average dimensions of the specimens along with 
the standard deviation are shown in Table 3-2. The adherend material properties used in Section 3.3.1 are 
representative of the adherends in this study, while the upper bound modulus is representative of the pure 
adhesive modulus. The specimens were placed in grips and pulled in tension under displacement control at a 
rate of 0.005 mm/s. Load and displacement data was recorded, and the joint strength was defined as the 
maximum load held by the joint. 

The grading was accomplished by means of inserting glass beads between the adhesive layers to change 
the properties of the adhesive. Bead placement was done purely by hand, and a photograph of the adhesive 
prior to doubler bonding of a representative specimen for each adhesive system is shown in Figure 3-15. 
Two specimens were made with no beads in the adhesive, two with a high density of beads uniformly 
spread throughout the adhesive region, and three with beads placed strategically by hand. Care was taken to 
ensure that the location of the stress concentration, where the two adherends meet, was devoid of beads with 
a gradual increase in beads going out from the center. 

An image of a single strap joint prior to complete joint separation is shown in Figure 3-16a. The 
resulting strengths of the joints tested are plotted in Figure 3-16b. Each bar represents a specimen; with two 
specimens with no beads, two with a uniform distribution of beads, and three specimens with graded beads 
in the adhesive. The specimens with no beads performed the worst, with an average strength of 3.18 KN. 
The uniform beaded joints had an average strength of 4.29 kN, and the graded bead specimens had an 
average strength of 5.43 kN. With so few specimens, no statistically significant conclusions can be drawn, 
but trends suggest that grading the beads could have resulted in higher strengths over the uniform 
specimens. 


Table 3-2. Geometric parameters for proof-of-concept testing. 


Iq (mm) I, (mm) t (mm) b (mm) 7 (mm) 
68.76+0.38 37.4940.25 1.08 +0.02 22.7340.17 0.42+0.03 
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Figure 3-15. Photographs of the adhesive layer before placement of the doubler for the single strap joint with 
(a) no beads, (b) uniform beads, and (c) graded beads. 
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Figure 3-16. Single strap joint test: (a) photograph of typical specimen prior to failure and 
(b) joint strengths of the specimens tested with different adhesive systems. 
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Figure 3-17. Post mortem photographs of the failed adhesive layer for the single 
strap joints with (a) no beads, (b) uniform beads, and (c) graded beads. 


A few words should be said about the effects of FGAs and failure modes. As was shown earlier in 
Figure 3-6, grading the adhesive has little or no effect on the adherends for the configuration studied. 
Therefore, a difference in strengths was only manifested when the failure occurred in the adhesive, not in 
the adherend or doubler. 

For these specimens, an improvement appears to exist using FGAs because the failures consisted mostly 
of cohesive failure with a few areas of adhesive failure (Figure 3-17). Another round of testing was 
performed with the same type of specimens cured at a lower temperature. These specimens all failed in the 
adherend, and there was no statistical difference in strengths between graded and non graded specimens. 
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Therefore, it should be noted that grading the adhesive can improve joint strength when failure occurs in or 
around the adhesive, and may even drive a joint to fail in the adherend when ungraded specimens fail in the 
adhesive, but FGAs most likely will not have a significant effect when failure already occurs in the 
adherends without grading. However, it is desirable for joints to fail in the adherend and adhesive at the 
same time so that the joint is maximized for strength, while minimized for weight. 


3.4 Concluding Remarks 


The peel stress concentration found in adhesively bonded joints can be reduced by grading the modulus 
of the adhesive material. More recently, grading of the joints in a continuous manner has become of interest. 
However, the actual use of such joints has been limited because of a few concerns which this study has 
addressed. These concerns were addressed by comparing the stresses in single adhesive, bi-adhesive, linear, 
and exponential functionally graded adhesives (FGAs) using an analytical model. This unique model allows 
efficient analysis of adhesively bonded joints in the context of a finite element framework without meshing 
the 2-D or 3-D details of the joint geometry. The efficiency of the model facilitated vast parametric and 
optimization studies, while still having the ability to be placed within a larger, global structural-scale model. 

The first concern was that gains to be made using FGAs were not worth the added manufacturing 
complications. It was shown that for the single strap joint configuration investigated, FGAs could reduce the 
maximum peel stress in the adhesive by up to 17% over the more compliant single adhesive joint without 
having an adverse effect on the adherend stress. Since the optimum grading for shear stress was not the 
same as that for peel stress, there was a slight increase in shear stress when peel stress was minimized. 
However, all grading functions resulted in a significant reduction (> 45%) in peel stress over the stiffer 
single adhesive joint without adversely affecting the load carrying capability or the stress in the adherends. 
Furthermore, decreasing the lower bound modulus of the grading caused even more dramatic stress 
reductions. Finally, it was shown that peel stress reductions comparable to those of the continuous functions 
studied here could be reached using a stepped function with around five or more different modulus 
adhesives, which is probably easier to manufacture. 

The second concern about using FGAs was that the flow of adhesive during manufacturing would 
change the shape of the grading and cause the grading to be ineffective. Therefore, a sensitivity study was 
conducted on the three FGAs to see the effect of changing the shape of the grading (I) and the lower bound 
modulus (E,) to reflect what might occur when adhesive is squeezed out of the joint. The maximum adhesive 
peel stress levels did increase for all FGAs when the lower bound modulus increased, but the exponential 
FGA proved to be quite tolerant to changes in grading shape. This could be the basis for a justification for 
using the more complicated exponential FGA over the bi-adhesive. Although the peel stress reductions were 
all very similar, the tolerability of the exponential grading to perturbation of grading shape which can occur 
due to adhesive flow might make it worth the extra complications. Also, it is recommended that care be 
taken to account for adhesive flow by adding more lower modulus adhesive at the end of the joint than is 
required, or using a grading system that does not change when the adhesive flows. 

Third, functionally graded materials used in structures can become useless or detrimental when the 
loading conditions are changed. It was shown that changing the loading configuration does not change the 
stress magnitude gradient of the adhesive in the case of a single strap joint. Consequently, the optimum 
grading shape and sensitivity was shown not to be significantly impacted by changing the loading 
conditions. However, in the case of moment and shear loading, both ends of the joint are subjected to high 
peel stresses, thus suggesting for generality that both ends be always graded. It was also shown that the 
stress magnitude gradient remains the same for other joint types; high stress at the adherend discontinuities 
with low stress elsewhere; thus suggesting the possibility of using the same adhesive grading for a myriad of 
joint configurations. However, additional parameters such as adherend stiffness, overlap length, and 
adherend thickness must be taken into account to conclusively determine the universality of a specified 
grading. 
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Additionally, proof-of-concept testing was conducted to show the potential of joints with FGAs. 
Grading was done by strategically placing glass beads to change the stiffness of the adhesive along the joint. 
Although not enough specimens were tested to result in significant conclusions, trends seem to indicate that 
these graded joints held over more load than specimens with no beads or uniformly distributed beads for 
failures dominated by cohesive and adhesive failures. Additionally, as expected, it was found that grading 
did not result in an increase in strength when the failure without grading occurred in the adherend. 

Through theoretical and experimental means, this study has shown the practical benefits of and offered 
some guidelines on how to effectively use and design joints with FGAs. FGAs have the potential to spread 
the stress more evenly across the adhesive without unrecoverable damage, making them another means of 
improving adhesively bonded joints. 
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Chapter 4 
Progressive Failure of Adhesively Bonded Composite Joints 


Enhanced finite elements are elements with an embedded analytical solution that can capture detailed 
local fields, enabling more efficient mesh independent finite element analysis. In Chapter 2, this method was 
applied to adhesively bonded joints. The adherends were modeled as composite Euler-Bernoulli beams, and 
the adhesive layer was modeled as a bed of linear shear and normal springs. The field equations were 
derived using the principle of minimum potential energy, and the resulting solutions for the displacement 
fields were used to generate shape functions and a stiffness matrix for a single joint finite element. In the 
current chapter, the capability to model non-linear adhesive constitutive behavior with large rotations is 
developed, and progressive failure of the adhesive is modeled by re-meshing the joint as the adhesive fails. 
The joint element is compared with experimental results for various joint configurations, including double 
cantilever beam and single lap joints. 


4.1 Introduction 


The joint element has been shown to predict adhesive stresses quite well under linear elastic conditions, 
but these conditions are not sufficient to predict joint strength. Modern polymeric adhesives are usually 
highly nonlinear, causing linear elastic analysis to be insufficient. Furthermore, the eccentricity of many 
joint configurations results in large rotations early on in the loading [44]. Large rotations are also necessary 
to predict the onset of imperfection-induced buckling, an aspect of great concern in thin-walled aerospace 
structures. 

Therefore, geometric nonlinear effects due to large rotations and material nonlinearity are pivotal in 
predicting the strength of a joint. This chapter will extend the previously created joint element to include 
these effects. Additionally, methods of growing an adhesive crack and adapting shape functions during the 
analysis will be presented in order to preserve the original intent of the joint element, which is to model a 
joint with very few elements. 

Currently, the scientific community seems to model the progressive failure of joints using fracture 
mechanics methods like LEFM and cohesive zone models, or continuum mechanics, both with dense-mesh 
finite elements [45]. Since the joint element is merely a tool, it will accommodate using inputs derived from 
either of these philosophies to govern the stress-strain relation of the adhesive. A method will be shown of 
characterizing the adhesive layer using either bulk adhesive tensile data as would someone using continuum 
mechanics damage progression, or fracture mechanics inputs like mode I strength and fracture toughness. 
The application of each will be demonstrated and results will be compared with experiments. 


4.2 Formulation 


The formulation of the joint element has been broken up into discrete parts, namely the co-rotational 
formulation, material nonlinearities, crack growth, adaptive shape functions, and adhesive constitutive 
modeling. Each section presents a formulation to address a certain aspect of the progressive failure of the 
joints. The co-rotational formulation addresses large rotations in joint problems while material nonlinearities 
show how nonlinear constituents are modeled. The crack growth formulation deals with the failure of the 
adhesive layer and adaptive shape functions are used to reduce the number of joint elements required to 
model nonlinear adhesive behavior. Finally, the last section illustrates a few methods of defining the 
properties of the adhesive based on several different experimental techniques. 


4.2.1 Co-rotational Formulation 


Consider a structure consisting of N layers of thin plates under cylindrical bending joined together by 
N-1 thin layers of a much more compliant adhesive material (see Figure 2-1). The plates are assumed to 
behave as “wide” Euler Bernoulli beams (hence the cylindrical bending assumption). The adhesive joining 
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the plates is modeled as a Winkler foundation. The plates can be isotropic, transversely isotropic, or a 
layered composite. The plates and adhesive are assumed to be under proportional loading, and are modelled 
as nonlinear elastic materials. 

A co-rotational formulation will be used to capture large rotations, and has been primarily adapted from 
prior work by Belutschko and Hsieh [46] and Crisfield and Moita [47]. This formulation tracks the rigid 
body rotation of an element through a local rotational coordinate system, and considers the rotations and 
deformations measured with respect to this rotated frame of reference to be small (Figure 4-1). The main 
benefit of this formulation is that the previously implemented code for the small rotation problem (see 
Section 2.2) can be utilized in subsequent calculations. 

The element has 2N nodes located at the boundaries of the centerline of the plates (numbered as shown 
in Figure 2-2), and the nodal displacements are defined as: 


F T 
q=[q' K q K q*] At 
where the superscript represents the adherend or plate, and 
qd=|a gadaa a fe 


refers to the horizontal, vertical, and rotational displacements of the left and right nodes in plate i 
respectively. 
4.2.1.1 Rigid Body Displacements 

The element has a local rotated coordinate system, X , which is rotated and translated relative to the 


fixed coordinate system, x , by an angle @ and a vector q, respectively (Figure 4-1). The translation and 


angle will be properly defined later. The nodal displacements of the element in the fixed coordinate system 
can be decomposed into rigid body displacements, q,,, , and displacements which only cause deformation 


in the body, q def ? in the relation: 


q es Grig + det : 4-3 


b) 


Figure 4-1. The nodal displacements can be broken up into two parts: a) rigid translation and rotations and 
b) local deformations. 
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The rigid body nodal displacements, q,,, , can be further decomposed into rigid body displacements 


resulting from rigid body rotation, q,., and displacements resulting from rigid body translation, q, : 


Giig =4: 749, 4-4 
The translational displacements, q,, are defined as: 
slat kK rr 
q, = [ a! qi | 4-5 


which is the horizontal and vertical displacements of the left node of the first plate and the rotation of the 
first adherend: 


a=[4 a 9) + 


Although the rotation is not necessarily part of the rigid body translation, it is more convenient to insert it 

into the translational rigid body displacements because each adherend will be rigidly rotated by the angle@. 
To find the rigid body displacements due to the rotation of the element about the first node, consider the 

right node of the ith adherend, node ir (Figure 4-2). Initially, node ir can be located relative to the first node 


by a position vector X,.. When the element rotates about the first node by the angle @, its new position 


relative to the first node can be expressed by an orthogonal transformation matrix as 1x where 


ir? 


c s O 
4-7 
T, =|-s 
0 
and s and c denote the sine and cosine of the angle @. Therefore, the displacement vector, qi, , of node ir 
due to rigid body rotation can be expressed as 
q" = (T’ —I)x, 
rot 1 ir. 4-8 


Figure 4-2. Displacements, qi, , of node ir due to a rigid body rotation of the joint element. 
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Translating this to all nodes and combining with Equation 4-4, the displacements due to rigid body rotation 
are 


Gig =q, + (T" rs I)x a 


where 


4-10 


and the vector X is simply a collection of the initial x and z coordiantes of the nodes, and is defined 
explicitly as 


X=[X, X,, K X;, X,, K Xvi ni va 


where the first subscript identifies the plate number, and the following letter, either / or r, refers to the left or 
right node respectively. The nodal coordinate vector for the ith adherend and the left node is defined as 


X,, =[X, Zi 0 sr 


while the coordinate vector of the right node is defined in an identical fashion. 


4.2.1.2 Determination of the Rotation Angle 


If the rotation is not constant within the joint, the rotation angle is an approximation. Adhering to the 
conventional approach for co-rotational beam formulations, the rotation angle was chosen to be the rotation 
of the first adherend as shown in Figure 4-1a. To find the transformation matrix of Equation 4-7, the sine 
and cosine of the rotation angle can be expressed as 


l 4-13 
s=sing=4 
L, 
and 
l 4-14 
c=cosg=— 
L, 
which are defined in terms of the nodal displacements of the first adherend by 
4-15 


11> 


Lise =1+q," qh 
L,=% -% 


where | is the original length of the element and 1, |,,, and I, refer to the current length of the 1* adherend 
and the length decomposed into x and z components Figure 4-3. 
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Figure 4-3. The initial and current lengths of the 1st adherend are used to determine the rotation angle. 


4.2.1.3 Local Coordinate System 


First, the internal force vector and stiffness matrix will be found in the local, rotating coordinate system. 
The stress and strain of the adherends and adhesive are assembled together in one stress and one strain 
vector as shown: 


A A 7m 7 in A i n WF 4-16 
6=[6, 6, K 6 Gy K Gy4 Gaya 6, | 


and 


i uk A. & a A” a 4-17 
e=[28, £4 MK € 8) K 894 Sey ey | 


where the overbar caret denotes quantities in the local rotating coordinate system, and the adhesive/adherend 
strain and stress vectors are as defined in Equations 2-2, 2-3, 2-12, and 2-13. Using beam theory and 
assuming small strains from the rotated coordinate system, the strains are related to the adherend centerline 
displacements, ti, by the equation 


4-18 


where the adherend centerline displacements are a collection of a vector of the centerline displacements of 
each adherend layer defined in Equation 2-33. Additionally, G is an assembly of the contributions of the 
adherend and adhesive layers assembled in the form 
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(s,] ~ 


[Gy] 


L Gay | 


[Gy] 


where the sub-matrices are defined in Sections 2.2.1 and 2.2.2. Furthermore, since the deflections in the 
local, rotated coordinate system are considered small, the shape functions derived for the linearly elastic 
case are used (Section 2.2.3). Using these shape functions, the local strain and displacements in the rotated 
coordinate system are related by the equation 


— Bq 4-20 
where B js defined as 
B=GN ror 


The principle of virtual work of the element can be written as: 


é(w'™ —-w*") - 0 4-22 
and the internal work can be written as the internal nodal forces multiplied by the nodal virtual 
displacements, or the integral of the strain energy density over the volume of the element: 

ow = oq't™ = [, w(q@av 4-23 


Where V is the volume of the element, and W(6q) is the strain energy density of the element resulting 


from a virtual displacement. Since the deformations are small relative to the rotated coordinate system, the 
internal virtual work can be rewritten as 


ow" =| dq'B'édVv as 


Assuming that external forces only occur as nodal forces and moments, the external virtual work of the 
element becomes 


ow = = oq't Ext 4-25 
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Finally, using Equation 4-22 and noting the fact that the virtual displacements are arbitrary, the resulting 
equilibrium equation is 


[ Bo av=t™ 128 
Now, the local internal nodal forces are 
f"o= | B'éaVv Pe 
V 
with the local stiffness matrix being given by 


k =| B’ DBdV 4-28 
Vv 


where D is an assembly of the adhesive and adherend matrices defined in Sections 2.2.1 and 2.2.2. Note 
that for linear elastic materials, the integration can be carried out analytically in the area, resulting in the 
relation given in Equation 2-47. 


4.2.1.4 Global Coordinate System 


Now we seek to find the residual and stiffness matrix in the global coordinate system. Since the internal 
work is not dependant on the frame of reference, one can write 


oq't Int _ oq't 4-29 


where the nodal virtual displacements in the global frame are related to those in the local rotated coordinate 
frame through the equation: 


oq = OG ig Be Toq 4-30 


making Equation 4-29 


oq f Int = (dq' dq 


T\pTg Int 4-31 
rig )T f 


Since rigid body motion does not result in the generation of internal forces, 


6q,,'T'f™ £(j 4-32 


and Equation 4-31 becomes 
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oq't Int _ oq Tt” . 4-33 


With the virtual displacements being arbitrary, the internal nodal force vector in the global coordinate 
system becomes 


gine = Tem 4-34 
To find the global tangent stiffness matrix, differentiation of Equation 4-34 gives 
5 Int _ oT fi +7! of 4-35 


The second term in the above equation becomes 


T" df" = T'kdq = T'kT6q-T'kT64q,,, | nee 


The last term vanishes because, as before, displacements resulting in rigid body translation and rotation do 
not generate any internal force. The first term on the right side of Equation 4-35 is somewhat more difficult 


to obtain. The difficulty lies in the fact that T contains sines and cosines of @, which in turn contain q 


and |. However, Crisfield [48] provides an approximation, which assumes that the extension Ol is small. 
Based on this assumption, the first term in Equation 4-35 can be rewritten as 


oT'¢™ = T,£'"0,5q 4-37 
where 
4-38 
Ts 
age O 
TT, 
and 
—s Cc 4-39 
Dga|se =s: 0 
0 0 O 
Similarly, 
= 4-40 
%, =|. 0K 0| 
and 
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1 4-41 
b4=7[s -c 0 -s c Of 


1 


Combining all of these equations, the global tangent stiffness matrix can be written as a combination of the 


material stiffness, K nat , and the geometric stiffness, k co in the relation 


k=k,,,, +k, 4-42 
where 
k,.. =T,'", 4-43 
and 
Kg = TRE a 


Both the geometric and material stiffness matrices are functions of the nodal displacements, making the 
system of equilibrium equations nonlinear. The Newton-Raphson method can be utilized to find the 
solution. It was already noted that one of the benefits of this method is that the formulation of the linear 
element can be utilized. Another major advantage of this method lies in the fact that the local rotational 
frame stiffness and internal force vectors are not functions of the nodal displacements. Since numerical 
integration occurs while finding these vectors/matrices, the integration must only be carried out once during 
the analysis. This saves a considerable amount of computational time, especially for an element like the joint 
element, which requires more refined integration for the higher order shape functions. 


4,2.2 Material Nonlinearities 


Since modern polymeric adhesives often display highly nonlinear material behavior, it was necessary to 
include material nonlinearities in the joint element to estimate joint strengths more correctly. A particularly 
simple nonlinear elastic stress law was chosen: 


6 =6(8) 4-45 


where the stress is some general function of the strain. The only major change from the previous co- 
rotational formulation is that Equation 4-28 becomes 


k =| B'D(@Bav a 


where the local stiffness matrix in the rotated coordinate system is now a function of the local 
displacements. 

Although it would be more correct to use an incremental flow type plasticity formulation that 
distinguishes loading and unloading stiffness, the simple nonlinear elastic relation, which assumes no 
permanent plastic strain, was chosen for several reasons. While this was chosen for simplicity sake, this 
decision can also be justified. The joint element is meant to be a design tool to give general approximations, 
so it is not expected that such a tool will be used in situations requiring unloading capabilities. Additionally, 
the nature of adhesively bonded joints is such that the high stresses occur in concentrated form at the joint 
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edges. Since the failing adhesive domain is eliminated in the iteration process (to be described later) the 
assumption of a nonlinear elastic type stress-strain law suffices for this modeling process since potential 
regions of “unloading” are minimal and contained in the regions which are eliminated. Thus, this 
assumption does lead to a meaningful rendition of the joint physics, yet facilitating an efficient (in the 
computational sense) solution strategy. 

One other aspect worthy of discussion is the integration requirements for the nonlinear material 
formulation. When the adherends have a nonlinear stress-strain relationship, Equations 4-46 and 4-27 must 
be integrated over X and Z at each Newton-Raphson iteration to allow a general stress-strain relationship. 
This causes a considerable increase in computational time. However, there are some cases when this is not 
necessary. If only the adhesive layers have a nonlinear stress-strain relation and adhesive model 1 (Section 
2.2.2.1) is used, integration over Z can be avoided because the stress is constant through the thickness of 
the adhesive layer. Additionally, if the functions for the nonlinear stress-strain relations are known (and 
simple enough), integration over Z can be accomplished analytically. However, this would mean that the 
formulation is only good for that specific stress-strain relation, and cannot be extended to other general 
relations. 


4.2.3 Crack Growth 


When some user defined failure criterion is reached in some part of the adhesive layer, that portion of 
the adhesive is considered “failed” and can carry no load and has no stiffness. Setting the stress and stiffness 
of that portion of the adhesive to zero is an easy way to model the failure of the adhesive, but the shape 
functions for the joint element were not originally calculated based on a joint with failed adhesive, and 
cannot accurately model this new situation. Therefore, as with more traditional shape function prescribed 
finite elements, more elements are required to accurately find the solution. In the case of failed adhesive, a 
great number of elements may be needed, as will be illustrated in Section 4.3.3. 

In order to increase the accuracy of the joint element after adhesive failure and crack growth, a method 
of removing the adhesive and adapting the mesh to the crack was devised. Since the joint element is meant 
to be used as a user defined element in a larger global assembly in commercially available finite element 
software, the mesh change would have to be strictly internal to the element so that the surrounding model 
does not have to change. Therefore, a sub-assembly method was devised to handle adhesive failure 
(Figure 4-4) and is outlined Figure 4-5. 


o inner nodes 
e outer nodes 


l crack 


(a) Uncracked (b) Partially Cracked (c) Fully Cracked 


Figure 4-4, Diagram showing a) an uncracked joint element, b) a partially cracked element, 
and c) a fully cracked joint element. 
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Figure 4-5. Flow chart showing how cracked element Sub-Assembly is incorporated 
into joint element solution procedure. 


First, when failure in the adhesive is detected, the element is replaced by a sub-assembly with three 
elements as shown in Figure 4-4. The length of the crack determines the lengths of the sub-assembly 
elements. Within a Newton-Raphson type solver, the nodal displacements are prescribed (guessed) and the 
stiffness and internal force vector for the element are calculated. These vectors/matrices for all of the 
elements in the assembly are assembled, boundary conditions and loads are applied, and the residual (error 
of the initial nodal displacement guess) is calculated. If the residual isn’t within some tolerated state, a new 
nodal displacement “guess” is calculated based on the previous displacement, residual, and stiffness values 
and the whole cycle repeats. 

In the case of a joint element with a crack, only the outer nodal displacements are prescribed since the 
global finite element assembly isn’t aware of the existence of the sub-assembly and the inner nodes. 
Therefore, the sub-assembly becomes a nonlinear model within another nonlinear model and must be solved 
with its own Newton-Raphson type solution procedure. The prescribed nodal displacements of the outer 
nodes become the boundary conditions for the sub-assembly, and the whole system is solved using a 
nonlinear solver. When the desired error tolerance is reached, a stiffness matrix and internal force vector for 
the sub-assembly has been calculated. However, these quantities still have the inner degrees of freedom 
contained within. The force vector and stiffness matrix are then reduced using the Guyan Reduction Method 
[49-51]. Once the internal degrees of freedom are removed, the stiffness matrix and force vector can be 
considered to be that of the equivalent joint element, and can be passed on to the global assembly. 

After the global system is solved, there is a check to see if the crack has grown, or if new adhesive 
failure has been detected. If this is the case, the sub-assembly is adjusted by changing the lengths of the sub- 
assembly elements, and the global system is re-solved. This is done until no new adhesive failure occurs and 
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the crack is in equilibrium. A crack scaling constant, C. , has been introduced to speed up or slow down 
crack growth as needed, and is used in the equation 


pour = you = Cc. d= a prev 4-47 


crack crack crack crack 


cur 
crac 


current crack length. Setting C. <0 causes the crack to grow further than detected, and is useful when 


where |?" is the previous crack length (prior to the global Newton-Raphson procedure) and | 


crack 


is the 


multiple iterations are needed to find crack equilibrium. Setting C, >0 causes the crack to grow less than 


detected, and is necessary when crack overshoot is a concern. 

The advantage of this method is that fewer elements are needed in order to accurately capture crack 
growth. One can use the minimum elements needed to accurately capture the material and geometric 
nonlinear effects without crack growth being a factor. This can mean dramatically reducing the elements 
required, especially when there is little material nonlinearity and strains in the joint are small. 

One of the major disadvantages of this method is the increased computational time. A local nonlinear 
problem must be solved within each iteration of the global nonlinear problem. Although the local nonlinear 
problem is always limited to three elements, it can significantly increase the runtime. Furthermore, the 
global load increment is repeated if the crack grows and the sub-assemblies need to be created or re-meshed. 
Although the crack scaling parameter can significantly help in limiting the iterations needed to find crack 
equilibrium this process can still be costly. However, the costs can be justified if joint strength prediction is 
of concern. Joint strength has been identified as a controlling factor in the ultimate load bearing capacity of 
many bonded structures. 


4.2.4 Adaptive Shape Functions 


When the adhesive in a joint element has a nonlinear stress-strain relationship as that shown in 
Figure 4-6a, the shape functions used earlier may no longer be well-suited because they were derived for the 
linear-elastic case. For highly nonlinear materials, many elements may be required to accurately capture the 
joint behavior. 

As the adhesive softens throughout the loading, the tangent modulus of the adhesive changes across the 
joint, as shown in Figure 4-6b. This resembles the functionally graded adhesive joints discussed in detail in 
Chapter 3. Using the non-constant coefficient ordinary differential solution procedure developed in Section 
3.2.1.1, the shape functions for the joint can be adapted to the adhesive softening at any point in the loading, 
which may decrease the number of elements needed. 
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Figure 4-6. For a single lap joint with (a) nonlinear material properties, the (b) tangent modulus of the 
adhesive layer becomes similar to a functionally graded adhesive in advanced stages of loading. 
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To this end, adaptive shape functions were applied to the joint element. After each load increment, the 
shape functions for the next increment are calculated for a joint with a functionally graded adhesive. The 
function for the adhesive modulus, F(x), is based on the tangent stiffness of the adhesive in the prior load 
increment. In this way, the shape functions adapt to the softening of the adhesive and can represent the 
behavior of the joint with very few elements (see Section 4.3.4). 

This has only been implemented for adhesive model 1 (Section 2.2.2.1) since the stress, and hence the 
tangent modulus, is a function of x only. However, a similar approach could be conceived for the other 
models by basing the adhesive modulus function on the adhesive centerline values (z,=0). Furthermore, a 
similar method could be applied to the adherends to improve elemental convergence when the adherends are 


highly nonlinear. 


4.2.5 Adhesive Model Characterization 


One of the most important inputs for determining the strength of a joint is the characterization of the 
adhesive constitutive response. There have been many methods of characterizing the adhesive material, but 
two have emerged as the most common: bulk adhesive tensile test and fracture mechanics characterization 
tests (DCB, ENF, etc.). Therefore, the following sections outline methods of using both bulk adhesive 
tensile test data and fracture mechanics inputs to characterize the joint. Ultimately, the test data available 
and personal preferences of the user will decide which route to take. 


4.2.5.1 Bulk Adhesive Tensile Characterization 

One common way of characterizing adhesive materials is by performing tensile tests on bulk adhesive 
specimens, such as those depicted in Figure 4-7. The following section will outline an approximate method 
for modeling the adhesive based on such adhesive characteristic data, and will discuss the formulation and 
underlying assumptions involved. 
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Figure 4-7. Adhesive may be characterized by (a) experimental bulk adhesive tensile tests, then 
(b) fitting a curve to the stress-strain plot. 
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Figure 4-8. Assuming that the adhesive is perfectly bonded to the adherends, 
the adhesive can be considered a constrained body under triaxial stress. 


If the adhesive is much deeper than thick (b,; << 7,;), it can be considered to be in a state of plane strain 
in the y direction (Figure 4-8), and the stress-strain relation shown in Equation 2-29 can be applied. 
Furthermore, if we assume that the adhesive is perfectly bonded to the adherends and that the adherends are 
much stiffer than the adhesive (E,; << E;), then it can be argued that the extensional strain in the adhesive is 
much smaller than the peel and shear components ( €,,, << &,,, V7, ) which is the root of the assumption: 


xai 


é.=0 4-48 


This assumption gives rise to the common practice in adhesive joint analysis of ignoring the extensional 
stress and strain in the formulation. With the extensional strains being relatively small, the strain energy of 
the adhesive layer is virtually unaffected and does not necessarily need to be included. Although the 
extensional strain is negligible, the extensional stress (in both the x- and y-directions) is not insignificant, 
placing the adhesive in a state of triaxial stress [52]. Equation 2-29 can be used to write the extensional 
stress in terms of the peel strain: 


6, =CVé 4-49 


Vv 4-50 


The same relation is true for the extensional stress in the y-direction, Oya . This can be used to find the 


extensional stress without necessarily including it into the formulation. Furthermore, the peel stress becomes 
a function of the peel strain only: 
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e Ly p 4-51 
0.3 = - ai€ai - 
(1 = 2V ai a + Vqi ) 


This shows that the effective “resistance” to deformation in the z-direction is amplified by a factor that 
depends on Poisson’s ratio. Although this relation is intended for linear elasticity, the relation was assumed 
to hold for the nonlinear stress-strain relation as well. Therefore, the stress-strain relation was redefined as: 


7 {= Vai 7 4-52 
Oi a ce ) 
(1 = 2V i a oF Vai ) 


which effectively increases the adhesive modulus. This relation can be used with adhesive models 1 and 2 
(Sections 2.2.2.1 and 2.2.2.2) to make them more accurate and to include the effects of having a triaxial 
stress state. 

A Von Mises failure criterion was chosen for this particular formulation, although the same formulation 
could easily be altered for a different criterion [53]. Applying the notation for the adhesive layer, assuming 
the shear stresses in the xy and yz planes to be negligible, and using Equation 4-50, the Von Mises 
equivalent stress in terms of the shear and peel stress components are 


4-53 
t=2V, 
: a Ou; +304, 
i> Vi 


vm 


Although the Von Mises equivalent stress is normally used to find the yield stress, in this case it will be 
assumed to hold for the entire nonlinear adhesive stress/strain response. Therefore, the Von Mises 
equivalent stress for a certain adhesive will be a nonlinear function of the adhesive strain found using bulk 
adhesive tensile tests (Figure 4-7): 


Vi 


Om = Faun = f (Epa) aaa 


To find the nonlinear curves approximating the peel and shear stress in the adherend, one more relation must 
be defined. A new variable will be introduced, y; , which represents the ratio of peel to shear stress for 


adhesive layer i of a particular joint configuration: 


allowing the shear stress to be defined as a function of the bulk stress: 
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2 4-56 


The method of finding the strain was a bit more arbitrary. Others have done this by utilizing a Von 
Mises strain criterion or similar methods [6,54,55]. For the current formulation, it was assumed that the bulk 
adhesive tensile specimen strain and the adhesive layer strains were linearly related to each other through 
the equations 


Voi = CoE wun =e 


and 


Eqi = CrE pu aap 
where the constants C, and C7 are found such that the initial slopes of the shear and peel stress-strain curves 
become the normal and shear modulus respectively. 


For an actual joint, the ratio of peel to shear stress, y; , not only varies across the joint, but changes 


during loading due to nonlinear geometric effects and nonlinear material effects. Therefore, this value will 
in actuality be a function of the joint geometry, loading, materials, and location within the adhesive in 
question. However, to simplify the determination of this value, it is proposed that one assume that the ratio 
of peel to shear doesn’t change significantly during the loading event and that only the stress at ends of the 
joint where the stress concentrations reside is important. The correctness of this first assumption will be 
tested later (Figure 4-33). Therefore, this value can be approximated by taking the ratio of the maximum 
peel to shear stress of the linearly elastic case as illustrated in Figure 4-9. 
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Figure 4-9. The peel to shear ratio for adhesive i can be approximated by dividing the 
maximum peel stress by the maximum shear stress for the linear elastic adhesive case. 
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For balanced joints with the same adherend materials and geometries, the maximum occurs on both 
ends of the adhesive and is identical on either end. However, for unbalanced joints, the stress concentrations 
at the ends of the adhesive can be of unequal magnitude. Finding the peel to shear ratio based on the higher 
and lower of the two stress concentrations can provide a good upper and lower bound to the nonlinear 
solution (see Section 4.3.4). 

In order to approximate the Von Mises failure criterion for uncoupled shear and peel, an uncoupled 
strain-based criterion was chosen that simply considered the adhesive failed when 


E aa Ve =1 4-59 
&, or 


where €. and /¥, are critical peel and shear strain values. These values are found by applying Equations 


4-57 and 4-58 to the maximum strain of the bulk adhesive tensile test data. 

Though it might seem unusual to use a strain-based criterion to approximate the Von Mises stress, it 
should be kept in mind that a Von Mises yield criterion was already applied to get from the bulk adhesive 
tensile test data to the peel and shear stress-strain relations. If the peel to shear ratio, Y, , was chosen 
correctly, both the shear and peel components should be close to their respective critical values at the same 
time. 


4.2.5.2 Fracture Mechanics Characterization 


The joint element model with adhesive model 1 is very similar to the Cohesive Zone Models (CZM) 
[44,56—58] and is inherently suited for fracture mechanics-type inputs. One of the main differences between 
most mainstream cohesive zone models and the joint element adhesive model lies in the thickness of the 
cohesive zone. Most CZM’s have no thickness, and lie at the interface between continuum elements. Since 
it has no thickness, a traction-separation law rather than a stress-strain law is defined for the CZM. Thus, 
cracks in the center of the adhesive layer can be differentiated from cracks at the interface by placing CZM 
elements at different locations within the adhesive, although this is computationally very costly. The joint 
element, on the other hand, resembles a cohesive zone with an explicit thickness. The entire adhesive layer 
is a single cohesive zone, and cracks in the middle of the adhesive are not differentiated from those at the 
interface. The traction-separation law can be transferred approximately to a stress-strain law by dividing the 
separation by the thickness as shown in Figure 4-10. 


Stress 


Strain 


Figure 4-10. Fracture mechanics properties such as critical stress 
and fracture toughness can be used to form an adhesive 
stress-strain law for the joint element. 
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For this type of adhesive characterization, the shear and peel responses are isolated and characterized in a 
series of experiments [59]. The peel and shear responses are considered to be uncoupled 

(see Equation 2-21), and depend solely on the vertical and horizontal separations of the adherends 
respectively (Equations 2-22 and 2-23). Typically, a critical stress and fracture toughness are identified for 
Mode I and Mode II. Since the joint element model does not have continuum elements to represent the 
adhesive, it is recommended that the initial slopes of the stress-strain laws be set to the elastic modulus for 
peel and shear. 


Finally, adhesive failure can be defined as occurring when 


i j 4-60 
G) (Se) 24 
Gic Gac 


The values of i and j can be chosen based on the preference of the element user. 
4.3 Results and Validation 


4.3.1 Geometric Nonlinearities 


To validate the co-rotational formulation, several example joint and beam configurations were analyzed 
using the joint element and compared with 2-D dense mesh finite element solutions with nonlinear 
geometric effects included or an analytical solution where applicable to demonstrate the joint element’s 
ability to capture large rotation situations and to demonstrate how many elements are typically required. 


4.3.1.1 Beam Problems 


The first example involves a class of beam elements, formulated like a joint element without an 
adhesive layer. The cantilevered beam is clamped at one end and displaced vertically at the other end in 
displacement-controlled loading. The beam is made of aluminum (E=70 GPa) and its shallow depth 
warranted a plane stress Euler-Bernoulli beam formulation for the joint element and 2-D plane stress 
elements for the dense mesh model. The geometric parameters are shown in Figure 4-11. A comparison of 
the load-displacement plots for the linear solution, 2-D dense mesh model, and joint element model is 
shown in Figure 4-12a. It can be seen that the joint element model agrees quite well with the 2-D dense 
mesh model. Additionally, a plot comparing the load-displacement curves for joint element models using 
different numbers of elements is shown in Figure 4-12b to show the number of elements required for a 
converged solution. Even with only two elements, the solution appears to be reasonably converged. 
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Figure 4-11. Geometric dimensions for aluminum beam problem for validation of the co-rotational formulation. 
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Figure 4-12. Comparison of (a) the linear solution, dense 2-D mesh, and joint element solution for the cantilever 
beam, along with (b) a convergence study for the joint element model. 
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Figure 4-13. Load vs displacement plots for a beam loaded in compression with a slight (a) moment and 
(b) imperfection angle. 


The next example features the same beam, except that it is compressed axially rather than loaded in 
bending to show how buckling loads can be found. Two different methods of applying imperfections are 
utilized. First, a slight moment is applied at the end of the beam (Figure 4-13a) and the load-displacement 
plots for different magnitudes of moments are compared with the well-known analytical solution [60]. 
Second, the beam is given a slight angle to induce buckling, and Figure 4-13b shows the results 
corresponding to different angles. In both examples, 50 joint elements were used and both agree well with 
the analytical solution, as expected. 


4.3.1.2 Single Adhesive Layer Joints 


The second type of example was single overlap joints, using the joint elements with a single adhesive 
layer and two adherends. These examples pointed more to the accuracy of the joint element formulation. 

The unbalanced single lap joint illustrated in Figure 4-14 was pulled in a displacement-controlled 
manner. The adherends were titanium (E=110 GPa) and aluminum (E=70 GPa), with EA 9394 as the 
adhesive layer (E=4 GPa, G=1.79 GPa). As before, the shallow width of the joint required the use of a plane 
stress joint element formulation and the use of 2-D plane stress elements for the dense 2-D finite element 
mesh model. The joint elements used adhesive model 1, and the Euler-Bernoulli adherend model. The joint 
element model had 40 beam elements with one joint element, while 154,000 elements were used for the 2-D 
dense mesh model (Figure 4-15a). A comparison of the load-displacement plots of the different models is 
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shown in Figure 4-15b. The joint element model was able to replicate the response quite well, even with 
only one beam element rather than 40. The actual joint region requires fewer elements because all of the 
bending takes place outside of the overlap region. The increased flexural rigidity of the overlap region 
causes it to rotate rigidly rather than bend. Therefore, more elements are required outside the overlap regions 
to capture the nonlinear geometric effects of the joint. 

The second single adhesive layer example is a layered beam subject to an axial compressive 
displacement. The adherends are aluminum, and the adhesive is EA 9394 (Figure 4-16). The width of the 
beam was 1 mm. Ten elements were found to produce a converged solution for the joint element model, and 
68,000 elements were used for the 2-D dense mesh model. The left end of the top adherend was clamped, 
while the left end of the bottom adherend was free to move vertically only. The right adherends were both 
constrained from rotating and a prescribed axial compressive displacement was applied. The beam was 
given an imperfection angle of 1° to induce buckling. 


e acat Titanium EA 9394 
aN "A es 


Aluminum 


Figure 4-14. Single lap joint used to validate joint element co-rotational formulation. 
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Figure 4-15. Comparison of (a) joint element and 2-D dense mesh finite element representation of the joint in 
Figure 4-14 and (b) the resulting load vs displacement plot. 
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Figure 4-16. (a) Geometric parameters and (b) buckled shape of layered beam under compressive loading. 
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Figure 4-17. Load vs. displacement plot comparing a dense, 2-D finite element 
mesh with the joint element using various adhesive and adherend models. 


As can be seen in Figure 4-17, the standard model with Euler-Bernoulli adherends and adhesive model 
1 buckled too early. To see whether adherend shear deformations played a significant role, Timoshenko 
adherends with adhesive model 1 was compared. Adherend shear effects played an insignificant role, 
making little difference. A model with Euler-Bernoulli adherends and model 2 adhesive was subsequently 
compared and showed a much better match with the 2-D dense mesh model. 

To explain why using adhesive model 2 over 1 made such a difference, one must remember the 
difference between the two models. Model 2 has the full shear strain to displacement relation, while the 
shear strain in model 1 is dependent on the difference in the difference in axial displacements of the 
adherends only. This is the difference between a Winkler Foundation with uncoupled shear and normal 
springs and a bed of shear and normal springs where the normal springs are attached to one another. This 
more accurate term for the shear shows how important shear effects are to buckling of sandwich-type 
structures. To further illustrate this point, the Young’s modulus and shear modulus of the adhesive were 
varied independently for the joint element model with Euler-Bernoulli adherends and adhesive model 1. The 
results of varying these parameters independently are shown in Figure 4-18. Increasing or decreasing the 
Young’s modulus by orders of magnitude has almost no effect on the response, while varying the shear 
modulus alone has a huge impact. Another notable observation is that there is a point where increasing the 
shear modulus doesn’t increase the buckling load. It can be assumed that this upper bound is the buckling 
load of the structure if it were a solid piece of aluminum. On the other side of the spectrum, when the shear 
modulus gets too low, the buckling load will approach two times the load of one adherend alone, as the 
adherends behave as independent beams. 
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Figure 4-18. Load displacement plots for the layered beam under compressive loading showing the effects of 
varying the adhesive (a) Young's modulus and (b) shear modulus illustrating the importance of the shear 
modulus on the buckling load. 
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Figure 4-19. Stress vs. strain approximation for aluminum beam depicted in Figure 4-11. 


4.3.2 Material Nonlinearities 


The material nonlinearity examples featured in this section only highlight adherend nonlinearity, rather 
than adhesive nonlinearity. Adhesive nonlinearity will be addressed in subsequent sections. This section 
only contains two brief examples, but the limitations and abilities of the joint element in modeling nonlinear 
adherends are shared by beam elements in general, and more in-depth discussion on these limitations and 
how to overcome them are dealt with extensively in literature [61-68]. 

The first example is the beam from Figure 4-11 where the adherend is modeled as being elastic- 
perfectly plastic with a yield stress of 300 MPa, as featured in Figure 4-19. Large rotations were not 
considered in this analysis to isolate the effects of the material nonlinearity. As can be seen from 
Figure 4-20a, the joint element compared quite well with the 2-D dense mesh finite element solution. 
However, as shown in Figure 4-20b, it took a good number of joint elements to converge to a solution. 
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Figure 4-20. Comparison of (a) the linear solution, dense 2-D mesh, and joint element solution for the cantilever 
beam with nonlinear materials, along with (b) a convergence study for the joint element model. 
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Figure 4-21. Plots showing the effect of element size on the load-displacement response of the single lap joint 
featured in Figure 4-14with nonlinear adherends for the (a) joint element, and (b) dense 2-d finite element mesh. 


The second example wasn’t quite so positive, and could serve as a good indication of why one would 
want to be careful when modeling joints with nonlinear adherends. The structure under question is the single 
lap joint shown in Figure 4-14, but with elastic-perfectly plastic adherends. The yield stress for the titanium 
was Set at 1050 MPa, and the aluminum was at 300 MPa. Large rotations were considered in the analysis, 
and the adhesive was given linear material properties to isolate the effect of nonlinear adherends. 

Figure 4-21a shows the load-displacement plot for the joint element model using different numbers of 
elements. As can be seen, none of the models are that far off of each other, but more elements are certainly 
necessary for a converged solution. However, the load-displacement plot did not resemble that of the 
Abaqus 2-D dense mesh model shown in Figure 4-21b. The load predicted by the 2-D dense mesh model 
drops after a peak, whereas the joint element model does not drop, but continues to hold more load. There 
are two explanations for this. First, since the stress-strain relation for the adherend is nonlinear elastic, 
unloading of the adherends is inaccurate. When the adherends first yield, the strain increases dramatically at 
one point (localization) while the rest of the adherend unloads. Since unloading is inaccurately captured in 
the joint element model, it continues to increase in load. The second discrepancy is that beam models still 
have the assumption that the displacement and strain vary linearly in the z-direction. Since this is not the 
case after yielding, the model is inaccurate after initial yielding. 


NASA/CR—2012-217606 73 


Another observation about the 2-D dense mesh finite element model is that the solution continues to 
change when the element size is reduced. This is due to the stress singularity at the reentrant corners. As the 
element size is decreased, the stress concentration rises and the adherends yield sooner and more 
dramatically. Furthermore, the solution cuts off after the peak for an element size of 0.15 and 0.1 mm. This 
cutoff was due to the commercial FE analysis software, which ends the analysis after the step size has 
become too small. This is also probably due to the stress singularity at the reentrant corners, and illustrates 
some of the potential difficulties of modeling joints. 

This last example illustrates why one should avoid using the joint element when failure of the joint is 
dominated by adherend yielding. It also brings out the need of applying some of the measures adopted for 
beam elements to the joint element to better capture the material softening of the adherends. 


4.3.3 Crack Growth 


To illustrate the benefits of growing a crack by re-meshing rather than just setting the failed adhesive 
stiffness and stress to zero, a bi-layered beam was pulled apart as shown in Figure 4-22a. The beam was 
5 mm wide, and the adherends had a stiffness of 100 GPa. The adhesive had a Young’s modulus of 1 GPa, 
and was linear up to failure, which occurred at 500 MPa. The simplistic linear-until-failure adhesive was 
chosen because an analytical solution can be found and because it allows crack growth without material 
nonlinearity, isolating this aspect of the joint element. 

Two different models were compared to show the benefits of re-meshing. First, rather than removing 
the adhesive and re-meshing, the stress and stiffness of the adhesive were simply set to zero when the stress 
reached 500 MPa. Second, the failed adhesive was removed and the element was replaced by a sub- 
assembly as illustrated in Figure 4-4b. The results of the two models with different ways of handling crack 
growth are shown in Figure 4-23. The benefits of re-meshing are clear. For the first model, the post-peak 
solution oscillates around the analytical solution with the oscillation amplitude reducing for more elements. 
The second model with the re-meshing, on the other hand, is extremely close to the analytical solution with 
just a single element. There is some oscillation after the peak, but this is suspected to be caused by crack 
overshoot. This effect, however, disappears entirely with only four elements. This example dramatically 
shows that re-meshing the element to represent crack growth can result in huge elemental savings over 
zeroing the adhesive stiffness. 
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Figure 4-22 Example of the peeling of a (a) layered beam where the adhesive is modeled as (b) linear until failure to 
demonstrate the joint element crack growth ability. 
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Figure 4-23. Load displacement plots for the peeling of a layered beam with different numbers of joint elements 
using (a) no re-meshing and (b) re-meshing. 


Table 4-1. For an unbalanced joint, the peel to stress ratio is different on each side 
and produces a different strength prediction. 


Side of Stress v Ou Tz Predicted 
Adhesive Concentration Oe) Mea) Strength (kN) 
Left Higher 1.63 29.0 Tel 710 
Right Lower 1.04 21.3 20.5 819 


4.3.4 Adaptive Shape Functions 


To show the effects of using adaptive shape functions, the single lap joint shown in Figure 4-14 with 
linear adherends and a nonlinear adhesive stress-strain relation is modeled with a dense 2-D mesh finite 
element model, the joint element model without adaptive shape functions, and the joint element model with 
adaptive shape functions. The adhesive had a bilinear bulk adhesive tensile test stress-strain relation similar 
to the plot in Figure 4-19 except that the linear properties are that of EA 9394 (E=4 GPa, G=1.5 GPa) and 
the bulk yield stress was 40 MPa. The procedure outlined in Section 4.2.5.1 was followed to find the peel 
and shear yield stress, O,,, and T,,,. The adhesive was allowed to yield indefinitely so that no crack would 


form or grow. This is an upper-bound prediction of joint strength according to the global yielding criterion 
proposed by Crocombe [69]. 

Since the joint was unbalanced, two peel to shear ratios were found; one on each side of the adhesive. 
The left side was the side with the greatest magnitude of adhesive stress, while the right side was a bit lower. 
Since the adhesive can yield indefinitely, the maximum load will not be reached until both sides of the 
adhesive begin to yield. Therefore, it was expected that the peel to shear ratio of the right side, the last side 
to yield, would result in the most realistic solution. The peel to shear ratios and peel and shear yield stresses, 
along with predicted joint strengths, are shown in Table 4-1. A comparison of the load-displacement 
response using the peel to shear ratio from the left (high ratio) and the right (low ratio) is shown in 
Figure 4-24a. This is expected to provide bounds for the solution. 
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Figure 4-24. Load displacement plots for the joint depicted in Figure 4-14 with an elastic perfectly plastic adhesive 
with yield stress of 40 MPa. Plot shows (a) results of basing the constitutive properties on the peel to shear ratio of 
the highest stressed side and the lower side, elemental convergence for the (b) 2-D dense mesh model, (c) joint 
element, and (d) joint element with adaptive shape functions. 


The load-displacement plot for different sizes of elements using the 2-D dense mesh model is shown in 
Figure 4-24b. As with the nonlinear adherend solution (Figure 4-21), the reentrant corners caused stress 
singularities, which cause the solution to be mesh dependent. However, it appears that for the element sizes 
shown, the joint element predictions provide an upper and lower bound for the 2-D dense mesh solution. 

The load-displacement results for the joint element model without adaptive shape functions are shown 
in Figure 4-24c. It takes a large number of elements to converge on a solution. In contrast, the joint element 
model with adaptive shape functions reaches a converged solution with just one element 
(Figure 4-24d), revealing the benefit of using adaptive shape functions. 


4.3.5 Experimental Comparison 


Two different joint types were modeled with the joint element to show some of its the capabilities and 
compare it to experiments. First, double cantilever beam (DCB) specimens tested by Song and Waas [70] 
were modeled with the joint element to validate the element with another finite element model and 
experimental data. The specimens were loaded and unloaded several times, which is the ideal situation to 
see the effects of assuming nonlinear elastic adhesive rather than elastic-plastic. Second, two additional 
DCB configurations were compared; one with a brittle adhesive and the other a much more ductile adhesive. 
Third, a single lap joint was modeled to show how the joint element compares to experiments in a mixed- 
mode test. The last two examples illustrate how bulk adhesive properties could be used to predict failure in a 
bonded joint. 
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4.3.5.1 Song and Waas DCB 


The ability of the joint element to predict the force vs. displacement behavior was assessed through 
comparison with experimental results published by Song and Waas [70]. The stable progressive failure 
exhibited in DCB specimens is ideal to show the capability of the joint element to fail progressively. Most 
single lap joints exhibit no stable crack growth; the crack grows almost instantaneously after crack initiation. 
Therefore, the DCB example illustrates the full capability of the joint element to grow a crack progressively. 
This particular data set was chosen because the authors stated that the failure was fully interlaminar and 
cohesive, which is the type of failure currently modeled by the joint element. 

Figure 4-25a shows a DCB specimen and the geometric parameters. The DCB specimens were 
constructed from two different 48 ply unidirectional composite laminates, E719/IM7 and E7T1/G40. The 
adhesive was one of the interlaminar matrix layers. The nonlinear stress-strain relation of the adhesive was 


based on G,, and E,,, shown in Figure 4-26, as was done by Song and Waas [70]. Normally, the mode one 


critical stress is preferred rather than the modulus, but this value was not provided by the authors. The 
geometric and material properties of the specimens are shown in Table 4-2. During the test, the loading was 
halted and the specimen was unloaded several times to measure the crack length within the specimen. The 
DCB specimen was modeled using one joint finite element accompanied by two beam elements on the top 
and bottom of the joint as shown in Figure 4-25b. 


Table 4-2. Material properties and geometric parameters for Song/Waas [70] DCB specimens. 


Adhesive Adhr. Geometric Parameters (mm) 
Ea Gic E 
Specimen (GPa) (N/m) (GPa) l a b t n 


E7T1/G40 4.1 335 116 200 52.6 15.5 4.65 35 
E719/IM7 3.3 1130 135 200 35.5 15.1 3.23 6 
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Figure 4-25. A typical DCB specimen, (a) the geometric parameters and boundary conditions for the DCB 
specimens and (b) the mesh for the joint element DCB model. 
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Figure 4-26. The stress-strain relation of the adhesive was defined based on the 
Mode I energy release rate and modulus. 
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Figure 4-27. Load vs. displacement curves for a) E719/IM7 and b) E7T1/G40 DCB specimens tested by 
Song and Waas [70] along with the present joint element model. 


Figure 4-27a compares the experimental force-displacement responses for the E719/IM7 DCB 
specimen acquired by Song and Waas [70], along with the response predicted by the present joint element 
model. As can be seen, the present model was quite accurate at predicting the progressive failure of the joint 
based on the given material properties and parameters. The behavior of the E7T1/G40 DCB specimens, 
shown in Figure 4-27b, was not predicted as accurately by the joint element model. The joint element model 
predicted a stiffer elastic response and higher peak load, but the subsequent response is captured quite well. 
This comparison shows that the joint element can be used with fracture properties similar to discrete 
cohesive zone models in predicting the behavior of DCB joints. 

Although these two DCB specimens were both loaded and unloaded several times, the non linear elastic 
material model used for the joint element adhesive was still sufficient to capture the overall behavior. This is 
because the stress concentration at the end of the joint causes the plastic zone in the adhesive to remain 
small. Therefore, the advantage of modeling the adhesive as elastic-plastic over non linear elastic is not 
great enough to justify the extra effort and complication. As long as the failure incurred due to crack growth 
is accounted for, the global response will be reasonably represented. 


NASA/CR—2012-217606 78 


4.3.5.2 Aluminum DCB Specimens 


DCB specimens were manufactured and tested to compare two paste adhesives, EA 9394 and EA 
9309.3NA. The results of these tests were used to assess the ability of the joint finite element to predict the 
difference in performance of two adhesives with very different stress-strain relations based on material 
properties obtained from tensile tests. Solid cylindrical specimens with a 3.175 mm diameter and 3.175 mm 
long test section were machined out of cast adhesive cylinders. The specimens were tested at NASA Glenn, 
and digital image correlation (DIC) techniques (Figure 4-28a) were utilized to obtain the axial strain of the 
specimen at different loads. The strains for several points in the gauge section were averaged, and the stress 
was found by assuming constant stress in the cross section. Figure 4-28b shows characteristic stress-strain 
data for the two adhesives, and the equations used to fit the data with a curve. The material parameters used 
in the equations are found in Table 4-3. The functions were chosen because they result in the same curve in 
compression and tension and seem to fit the data adequately. 


Table 4-3. Material properties and geometric parameters of DCB specimens. 


Adhesive Adhr. Geometric Parameters (mm) 
Eq Ou E 
Specimen (GPa) (MPa) — &qaii v (GPa) l a b t n 
EA 9394 4.2 496 0.016 0.4 69 152.4 63.5 25.4 12.7 0.6 
EA 9309.3 
NA 2.7 41.3 0.068 0.42 69 152.4 63.5 25.4 12.7 0.55 
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Figure 4-28. Using a) cylindrical tensile specimens and digital image correlation, the b) stress-strain relations of the 
adhesives EA 9394 and EA 9309.3NA could be defined by curve-fitting experimental data. 
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Adherends Adhesive 
Figure 4-29. Two adherends of an EA 9394 DCB specimen after complete failure. Adhesive found on both 
adherends indicates that failure occurred within the adhesive layer as desired. 


Since the joint element considers only failure in the adhesive (cohesive failure), care was taken to 
ensure that the interface between the adherends and adhesive of the DCB specimens would not fail. The 
adherends were 7071 T6 Aluminum, and the surfaces to be bonded were sanded, etched in lye, and anodized 
in a sulfuric acid solution prior to bonding [71]. This treatment was sufficient to produce failures in the 
adhesive layer, as can be seen on the failure surfaces of a post-mortem specimen in 
Figure 4-29. The failed specimen has adhesive covering both adherends, which means that the interface was 
not the plane of failure. Glass beads were used to maintain a consistent bond line thickness throughout the 
specimen, and pressure was applied to the specimen during curing. The specimens were allowed to cure for 
seven days at room temperature. 

Three DCB specimens for each adhesive were tested on an Instron machine at 0.5 mm/min. All 
specimens failed cohesively like the specimen in Figure 4-29. The load-displacement curves for all six 
specimens are shown in Figure 4-30. The high strain-to-failure of EA 9309.3NA caused these specimens to 
hold over two times load of EA 9394 specimens. The EA 9394 specimens exhibited a load plateau rather 
than dropping in load after adhesive failure was initiated. It is possible that air bubbles in the adhesive 
caused the adhesive to fail prematurely, allowing the joint to not drop in load carrying capacity after failure 
initiated. 

The joint element model was able to capture the behavior of the joints rather well. It was found that 
compliance in the experimental load train caused the models to over-predict even the initial linear portion of 
the loading. To compensate for this system compliance, the length of the adherends was increased by 7.5% 
for both DCB specimen types. This number was determined by fitting a linearly elastic model to the initial 
portion of the experimental force/displacement plot. These experiments were very effective in displaying the 
ability of the joint element to predict failure, along with showing how constitutive relations can be applied to 
get progressive failure of a thin adhesive layer. 


4.3.5.3 Harris and Adams Single Lap Joint 


Lastly, the joint element was compared with experimental data published by Harris and Adams [55] on 
single lap joints. The tests were carried out according to ASTM D1002-72 specifications. The geometric 
parameters are shown in Figure 4-31. The adhesive was MY750 and three different aluminum alloys served 
as the adherends. The only difference between the alloys was the 0.2% proof stress, as shown in 

Table 4-4. The adherends were modeled with an elastic-perfectly plastic stress-strain relation similar to 
the aluminum response shown in Figure 4-19. The adhesive, MY750, was characterized using bulk adhesive 
tensile tests, and the bulk adhesive stress-strain relation is shown in Figure 4-32a. 
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Figure 4-30. Load vs. displacement curves for DCB specimens with aluminum adherends and EA 9394 and EA 
9309.3NA adhesive, along with the joint finite element model prediction. 
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Figure 4-31. Geometric parameters for single lap joint tested by Harris and Adams [55]. 


Table 4-4. Material properties of the single lap joint adherends and adhesive [55]. 


E (GPa) v 0.2% Proof Stress (MPa) 
MY750 3.44 0.4 - 
Aluminum 2L73 70 0.34 430 
Aluminum BB2hh 70 0.34 220 
Aluminum BB2s 70 0.34 110 
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Figure 4-32. (a) Stress-strain relation for bulk adhesive, along with peel and shear components for a single 
lap joint with y=1.4, and (b) corresponding load-displacement plots. 


Table 4-5. Experimental and predicted strengths of the single lap joint. 


Adherend Experimental Strength (kN) Predicted Strength (kN) 
2L73 4.8 + 0.57 4.46 
BB2hh 5.0 + 0.38 4.52 
BB2s 3.5 + 0.32 5.00 


The method outlined in Section 4.2.5.1 was followed to find the adhesive peel and shear stress-strain 
relation. First, the joint was analyzed with linear material properties and small rotations, and the peel to 
shear ratio, Y was found to be 1.4. Using this value, the Young’s modulus, and the Poisson’s ratio, the bulk 


adhesive tensile data was converted to the peel and shear stress-strain relations shown in 

Figure 4-32a. Using this, the joint was modeled with 20 beam elements and one joint element and was 
loaded in a displacement controlled manner until the peak load had been reached. The load-displacement 
plots for the single lap joints with different aluminum alloys are shown in Figure 4-32b, and the results are 
compared with the experimental values found by Harris and Adams [55] in Table 4-5. 

The joint with the 2L73 adherends failed without the adherends reaching the yield stress, while the 
BB2hh adherend joint had small amounts of adherend yielding and the BB2s joint was dominated by the 
effects of adherend yielding. Looking back to the single lap joint example of Section 4.3.2, the point was 
made that adherend yielding is not accurately captured by the current formulation of the joint element. As 
expected, the specimen with no signs of adherend yielding, 2L73, had a predicted strength well within the 
experimental error. The specimen with slight yielding, BB2hh, had a predicted strength slightly outside of 
the error range of the experiment. Finally, the BB2s adherend joint, being totally dominated by adherend 
yielding, had a predicted strength much higher than the experimental value. However, if one again uses the 
single lap joint of Section 4.3.2 as an example, one could easily imagine that if adherend plasticity were 
accounted for in a more accurate manner, the predicted peak load would be somewhere around the elbow 
where the slope first drops, around 3 kN. This would bring the prediction much closer to the experimental 
value. Unfortunately, as predicted in Section 4.3.2, the joints with more adherend yielding predict strengths 
increasingly deviating from the experimental value. 

If the elbow is taken to be the point of failure for the BB2s specimens, all three predictions would be 
lower than the experimental strength. There are several possibilities for this discrepancy. The first is that the 
actual joints had quite sizeable fillets at the ends of the adhesive. Although it has been shown that spring- 
type joint models, like the joint element, predict stresses within the bondline similar to those in joints with 
fillets [6], the fillet might reduce the stress enough to increase the strength slightly. Furthermore, the peel to 
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Figure 4-33. Peel to shear stress ratio in adhesive layer of the single lap joint as a function of end displacement. 


shear ratio, y, was only approximated base on the linear elastic joint. However, large rotations and the 
accompanying nonlinearities change the peel to shear ratio, making it a function of the loading. Figure 4-33 
shows the value of y as a function of the end displacement, 4. It can be seen that the peel to shear ratio 
drops early on in the loading. Therefore, y could be adjusted to yield a more accurate answer. 

This comparison showed that, as expected, the joint element is less than accurate with regards to 
adherend material nonlinearity. However the method devised to use bulk adhesive tensile data appears to 
have been successful in approximating the strength of this single lap joint. For most advanced composite 
joints, the adherends display brittle failure, so capturing adherend yielding is of secondary importance. 
However, a more precise model could be implemented to consider adherend damage. 
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Chapter 5 
Conclusion 


In conclusion, a brief summary will be provided for the dissertation and the main points and conclusions 
will be highlighted. Following the summary, recommendations for future work to improve and expand the 
current models will be suggested. 


5.1 Brief Summary 


In Chapter 2, a formulation was laid out to find the stiffness of a joint element under linearly elastic 
conditions. The joint element was created by finding the stiffness of a beam-on-elastic foundation type 
model, with any number of beams stacked and attached together by adhesive layers. The adherends were 
considered beams under cylindrical bending, and their displacement fields were reduced down to a function 
of x only. The adhesive was assumed to be linear in the transverse direction and written in terms of the 
adherend displacements. Using such a structural model allowed the governing equations to take the form of 
a system of ODEs, which were subsequently solved in a semi-numerical fashion to yield the stiffness 
matrix. The formulation was done in a general fashion, which allowed several different adherend and 
adhesive models to be inserted. Additionally, since modern composite joints are often very complex with 
tapers and other features, a basic “building block” framework was created to allow the modeling of realistic 
joints. 

Several different joint configurations were used to validate the joint element by comparing them with 2- 
D dense mesh finite element models. First, a parametric study compared several of the adhesive models for 
different joint geometries. It was found that the simplest model, which was primarily used for subsequent 
analysis, provided good results for thin, long adhesive layers. Additionally, more complex joints with steps 
and tapers were used for validation, illustrating the ability to model such joints. Finally, a joint with 
laminated composite adherends was modeled and found to agree well with a model generated by 
commercial sizing software. 

Chapter 3 veered away from the formulation of the joint element, and demonstrated its usefulness in 
design and iteration type analysis. The model was used to show the practicality of using functionally graded 
adhesives in advanced composite joints to spread out the stress concentration and increase the overall joint 
strength. This was done by addressing concerns expressed in the adhesive bonding community, namely 
whether the gains are worth the increased complication, whether the flow of adhesive during manufacturing 
can neutralize or even negate the benefits of grading, and whether the same grading can be used for different 
loading scenarios. 

These concerns were addressed by modeling a sample joint configuration with different types of grading 
functions. The analysis showed that adhesive stress reductions up to 17% were possible for the 
configuration studied. Furthermore, a sensitivity study was conducted to show the effect of changing the 
grade shape, which simulated the flow of adhesive during curing. It was shown that for some types of 
grading, the change in shape had only a minimal effect on the beneficial aspects of the grade. However, this 
was not true for all grading functions. Last of all, research revealed that for the configuration studied, the 
optimal grading did not change by changing the loading, nor did the optimum change for different joint 
types. This finding suggests that the same gradation could be used universally for any number of joint types, 
making a case for the mass-production of functionally graded adhesive tapes. All-in-all, the joint element 
was very useful for exploration of functionally graded adhesive joints in an extremely efficient manner. 

As an additional demonstration, joints with functionally graded adhesives were manufactured and tested 
to failure. Although the methods of grading were not precise, this demonstrated the potential for such joints. 
With test results alone, it was found that the graded adhesive specimens had an average strength increase of 
10% over uniform adhesive specimens. 
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Finally, in Chapter 4, the linear elastic joint element concept was extended to include large rotations, 
material nonlinearity, and adhesive failure. Large rotations, which occur commonly in adhesively bonded 
joints, were handled through a co-rotational formulation. This formulation separated the displacements into 
rigid body displacement and local deformations about some rotated local coordinate system. The local 
deformations are assumed to be small, so the linear formulation can still be used. Material nonlinearities 
were included into the formulation. However, a nonlinear-elastic model was adopted for simplicity. It was 
shown through examples that, while this model was sufficient for the adhesive layers with high stress 
concentrations and often small plastic zones, it was not accurate for a description of the adherend materials, 
especially in the post-yielded stated. Problems arise with excessive adherend yielding and it is suggested the 
modeling of such joints with the joint element be avoided. On a positive note, the joint load associated with 
adherend yielding can be viewed as an upper limit load for the structural joint, predicted using the joint 
element. 

Adhesive failure and crack formation and growth were accounted for through an internal re-meshing 
process. The element with an internal crack was replaced by a sub-assembly with the failed adhesive 
removed. This method added to the computational steps that needed to be taken during the analysis, but 
decreased the number of elements needed to capture progressive failure. To further decrease the number of 
elements required, a method of adaptive shape functions was employed. This used the formulation for the 
functionally graded adhesive element to solve for new shape functions at each load step based on the tangent 
modulus of the adhesive in the prior load step. For highly nonlinear adhesive materials, this method results 
in a dramatic reduction of elements needed. 

Finally, methods of finding the nonlinear peel and shear stress-strain curves for the adhesive based on 
experimental procedures were outlined. First, using bulk adhesive tensile data, the response was broken up 
into shear and peel components for a certain joint configuration. This allowed the adhesive to be 
characterized with one test, but limited the shear and peel characterization to be specific to a certain joint 
type, geometry, and materials. Next, the resemblance of the adhesive model to cohesive zone models made 
it a natural candidate for fracture properties such as strength and fracture toughness. Tests are conducted to 
isolate the shear and peel “modes” and characterize them separately. This has the disadvantage of requiring 
more tests, but seems to have fewer assumptions involved. The joint element was constructed such that the 
user can use whichever way he/she feels is the most correct. 

To cap off the study, the joint element was compared with several experimental data sets, including 
double cantilever beam and single lap joints. It was found to compare well when little or no adherend 
yielding was involved. The joint element was shown to have practical use in estimating the strength of an 
adhesively bonded joint. 


5.2 Future Work 


As with any worthy research topic, there is a vast amount of possibilities for further work to either 
extend the concepts or models outlined here or improve upon their current state. Here is a list of some of the 
main areas for future work. 


5.2.1 3-D Plate Joint Element 


The current joint element formulation used a beam type formulation, assuming plane stress or plane 
strain through the thickness. To make the joint element really useful for vehicle and structural designers, a 
3-D plate or shell type element needs to be constructed. This would require solving a system of PDE’s 
instead of ODE’s, and some way would need to be devised to convert the nodal displacements into edge 
boundary conditions. If this were done, the joint element would be truly ready for distribution and 
widespread use in the engineering community, and its utility would be many-fold. 
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5.2.2 Joint Element Extensions 


The general formulation of the joint element allows it to be extended into other applications. One of the 
easiest and most natural extensions is to create a Sandwich Element that can be used to analyze sandwich 
structures. Similarly, with the ability of having any number of “adherend” layers stacked on top of each 
other, the element could be used to find interlaminar stresses in a composite laminate, and even model 
delamination. If the 3-D plate joint element were created, a composite shell element would be extremely 
useful for designers in predicting the delamination of thin-walled composite structures. 


5.2.3 Functionally Graded Adhesives 


The study presented here on the usefulness of functionally graded adhesives was only an introduction 
and viability study to show that this concept warrants attention. Further work in this area would involve 
identifying a practical grading method which has the ability of providing very exact higher order gradations 
and demonstrating even greater benefits than the preliminary specimens shown herein. Current advances in 
nanotechnology may provide insights into manufacturing, in an automated manner, functionally graded 
polymer adhesives [72]. Furthermore, the adhesive would need to be characterized fully for the different 
grades. Then, the joint element model could be in conjunction with optimization methods to design the most 
beneficial grading functions. If the concept proves successful, mass-production would be the next logical 
step. 


5.2.4 Adherend Yielding 


Finally, the current joint element model showed a deficiency in its ability to represent adherend 
yielding. If joint performance beyond adherend yielding is required, then the adherend constitutive model 
would need to be represented through an appropriate incremental flow theory of plasticity approach with 
softening hinges [63], however, at the expense of computational efficiency. 
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